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Abstract 


Modern  fighter  designs  have  been  associated  with  lateral  self-excited  oscillations  known  as 
“Wing  Rock”.  Simulations  of  wing  rock  by  flow  visualization  utilities  have  been  encouraged  to 
develop  a  complete  understanding  of  the  fluid  mechanism  that  drives  the  motion  as  well  as  for 
prediction  purposes.  Previous  wind/water  tunnel  simulations  have  been  limited  to  a  single  degree- 
of-freedom  (DoF)  in  roll.  Numerical  simulations  utilizing  computational  fluid  dynamics  (CFD) 
have  also  been  limited  to  a  single  DoF  in  roll  and  applied  to  simple  flying  delta-wing 
configurations.  The  loss  of  simulation  accuracy  due  to  reducing  the  actual  wing  rock  six  DoF  to  a 
single  roll-only  DoF  has  not  as  yet  been  fully  investigated.  In  this  study,  for  the  first  time,  wing 
rock  is  computationally  simulated  in  three  DoF:  roll,  sideslip,  and  vertical  motion  to  study  the 
effect  of  adding  the  sideslip  and  vertical  motion.  The  results  are  for  a  generic  fighter  model 
consisting  of  a  fore-body,  a  cropped  delta  wing,  and  a  vertical  fin.  The  effect  of  including  the 
vertical  fin  is  also  studied. 

Computational  simulation  of  wing  rock  in  three  DoF  for  a  full  generic  fighter  model  is  a 
multidisciplinary  problem  involving  grid  generation  (which  requires  the  utilization  of  the  multi¬ 
overset  mesh  technique),  aerodynamics,  and  rigid-body  dynamics.  The  coupling  of  the  flow 
governing  equations  and  the  rigid-body  dynamics  equations  has  been  done  utilizing  an  innovative 
sub-iterative  scheme  that  works  on  the  coupled  set  of  equations.  The  advantage  of  implementing 
this  sub-iterative  scheme  is  to  reduce  the  inherent  time  lag  between  the  solutions  of  the 
aerodynamics  and  rigid-body  dynamics. 

The  interaction  of  aerodynamics  and  rigid-body  dynamics  during  a  single  DoF  wing  rock  for 
the  wing-body  configuration  has  been  studied  via  snap  shots  of  a  cross-plane  stagnation  pressure 
distribution  and  tracing  the  instantaneous  locations  of  vortex  burst  for  an  entire  cycle  of  wing 


xv 


rock.  An  innovative  explanation  of  the  fluid  mechanism  that  drives  and  sustains  the  motion  has 
been  introduced. 

The  effect  of  adding  the  sideslip  and  vertical  motion  DoF  to  the  simulations  of  the  wing- 
body  configuration  was  found  to  delay  the  onset  and  to  reduce  the  amplitude  of  wing  rock  by 
about  50%  with  surprisingly  no  change  in  frequency.  In  an  effort  to  correlate  the  variation  of  the 
rolling  moment  coefficient  with  the  associated  sideslip  for  a  full  cycle  of  wing  rock,  a  pseudo 
dihedral  derivative  was  found  to  have  a  clearly  stable  mean  value  (-0.0591  rad"1). 

The  wing  rock  simulation  in  three  DoF  was  repeated  for  the  full  generic  fighter  model  with 
the  fin  included.  The  aerodynamic  effect  of  the  fin  was  found  to  significantly  delay  the  vortex 
burst  on  the  upper  surface  of  the  wing.  The  net  effect  of  the  fin  was  found  to  augment  the 
damping  of  the  oscillations  with  significant  increase  in  frequency. 
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SIMULATION  AND  ANALYSIS  OF  WING  ROCK  PHYSICS 
FOR  A  GENERIC  FIGHTER  MODEL  WITH 
THREE  DEGREES-OF-FREEDOM 


Chapter  1-  Introduction  and  Background 


1.1  Introduction 

Present-day  fighters  are  tailored  specifically  to  operate  at  high  speeds  and  high  angles  of 
attack  to  meet  the  demands  for  air  superiority.  A  configuration  consisting  of  thin,  low  aspect 
ratio,  highly  swept-back  or  delta  wing  and  a  long  slender  fuselage  has  been  adopted  in  most 
modern  fighters.  Also,  aerodynamic  features  and  control  systems  have  been  devised  to  provide  a 
significant  increase  in  the  operational  angle  of  attack  range.  The  employing  of  these  modern 
configurations  led  to  a  high  reduction  in  drag,  especially  in  the  transonic  and  supersonic  regions, 
and  consequently  a  considerable  increase  of  the  maximum  attainable  speeds.  The  concentration  of 
mass  around  the  longitudinal  axis  has  significantly  reduced  the  rolling-axis  moment  of  inertia, 
and  thus  augmented  the  maximum  attainable  roll  rates.  Two  key-elements  of  air  superiority, 
speed  and  maneuverability,  have  been  dramatically  improved  in  these  modern  configurations. 
However,  on  the  downside,  the  low-speed  aerodynamic  characteristics  have  generally 
deteriorated.  Flying  these  configurations  at  low  speeds  requires  relatively  high  values  of  angle  of 
attack,  which  promotes  a  typical  three-dimensional  vortical  flow.  Highly  swept-back  or  delta 
wing  configurations  are  known  to  have  stronger  directional  stability  compared  with  its  lateral  one. 
This  reduction  in  lateral  stability  together  with  the  reduction  in  the  rolling  inertia  increases  the 
susceptibility  of  modern  configurations  to  rock  under  any  lateral  perturbations.  Once  a  lateral 
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perturbation  is  triggered,  the  vortical  flow  starts  to  interact  with  the  body  leading  to  oscillations, 
primarily  in  roll,  that  vary  in  nature  according  to  the  value  of  angle  of  attack.  At  low  angles  of 
attack  the  resulting  oscillations  are  damped,  but  at  higher  angles  of  attack  the  damping  is  reduced 
and  a  limit  cycle  may  be  reached.  The  limit  cycle  condition  is  commonly  referred  to  as  “Wing 
Rock”,  which  has  been  investigated  through  intensive  research  in  the  past  two  decades. 

Lateral  position  precision  is  of  prime  importance  for  handling  quality  assurance  of  a 
modern  fighter,  and  thus  it’s  rating.  Target  tracking  (Roll-and-Hold),  air  refueling  (receiver),  and 
landing  in  a  crosswind  condition  are  examples  of  maneuvers  where  stable  lateral  position  control 
is  vital.  For  example,  in  the  takeoff/landing  of  the  F-14A,  at  28-deg  angle  of  attack,  the  stall  is 
characterized  by  divergent  wing  rock  that  can  reach  90  degrees  within  6  seconds  if  the  stick  is 
held  back  [2:876].  Moreover,  high-frequency  wing  rock  can  lead  to  a  phenomenon  known  as 
“Pilot-Induced  Oscillations”  or  PIO,  if  the  pilot’s  response  to  control  (or  suppress)  the  oscillations 
is  slower  than  the  motion. 

One  of  the  research  interests  in  the  area  of  wing  rock  is  simulation,  which  provides  useful 
information  in  the  preliminary  design  phases  about  limit  cycle  characteristics  and  the  onset  angle 
of  attack.  Moreover,  simulation  can  be  used  in  the  analysis  and  study  of  the  physical  mechanism 
that  drives  the  motion.  Simulations  of  wing  rock  began  experimentally  in  wind/water  tunnels  for 
a  simple  delta  wing  [24,31,43]  with  a  single  degree-of-freedom  (DoF)  in  roll  and  later  were 
extended  to  full-fighter  models  with  a  single  DoF  [4,42].  With  the  computer  revolution  of  the 
80s,  numerical  simulations  of  wing  rock  have  been  encouraged.  Computational  simulations  have 
been  shown  to  provide  results  that  agree  well  with  experimental  data.  The  developments  in  the 
computational  resources  have  led  to  considerable  improvements  in  the  accuracy  of  the 
computational  simulations.  Generally  though,  computational  simulations  have  been  limited  to 
delta-wing  configurations  with  a  single  DoF  in  roll. 

The  aircraft  degrees-of-freedom  are  in  general  coupled;  where  the  degree  of  coupling 
depends  on  the  nature  of  the  maneuver  and  flight  conditions.  One  of  the  most  useful  studies  in 
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flight  dynamics  is  to  determine  for  a  particular  maneuver,  either  theoretically  or  experimentally, 
the  dominant  degrees-of-freedom,  or  equivalently  a  reduced  order  model  (ROM).  The  use  of  a 
ROM  provides  a  very  close  response  to  the  original  model  at  lower  computational  or 
experimental  efforts.  In  complex  maneuvers,  it  is  hard  or  even  impossible  to  theoretically  prove 
that  the  ROM  is  adequately  representing  the  original  model.  So,  for  most  of  the  complex 
maneuvers,  experiments  are  used  to  prove  the  reliability  of  the  ROM.  In  the  latest  reviews  on 
wing  rock,  Katz  [27]  and  Mabey  [37]  have  reported  (results  of  which  will  be  presented  in  detail  in 
Chapter-2  that  previous  simulations  of  wing  rock  have  been  conducted  with  only  a  single  DoF  in 
roll.  Also,  the  delta  wing  has  only  been  modeled  to  represent  the  aerodynamics  of  the  full  aircraft 
configuration.  The  reliability  of  the  computational  and  experimental  simulations  based  on  this 
ROM  has  not  as  yet  been  fully  investigated. 

In  previous  work  by  Liebst  [33,34,35]  on  the  F-15,  the  sideslip  has  been  shown  as  the 
second  major  contributing  DoF  after  roll.  In  a  study  by  Saad  and  Liebst  [48]  on  the  F-15,  the 
dihedral  derivative,  and  hence  the  sideslip,  has  been  shown  as  the  major  contributing  derivative 
that  controls  the  value  of  wing  rock  onset. 

In  the  present  study,  the  accuracy  of  the  prediction  of  the  onset  angle  of  attack  and  the 
limit  cycle  characteristics  based  on  a  simple  delta-wing  with  a  single  DoF  in  roll  will  be  analyzed. 
The  simulations  of  the  wing  rock  oscillations  for  a  generic  wing-body  configuration  modeled 
with  a  single  DoF  in  roll  will  be  compared  and  analyzed  with  that  corresponding  to  three  DoF 
modeled  with  roll,  sideslip,  and  vertical  motion.  Also  the  simulations  will  be  done  for  the  wing- 
body  with  and  without  the  vertical  fin  to  study  the  effect  of  including  the  fin  in  the  prediction 
studies. 
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1.2  Various  Types  of  Wing  Rock 

Vortical  flow  over  highly  swept-back  or  delta  wings  at  high  angles  of  attack  is  highly 
nonlinear  due  to  vortex  asymmetry,  streamwise  and  cross-flow  separation,  and/or  vortex 
breakdown.  Due  to  the  complexity  of  the  interaction  of  one  or  more  of  these  phenomena  with  the 
rigid-body  motion,  many  studies  have  been  conducted  to  thoroughly  understand  the  physics  of  the 
phenomena. 

Previous  research  [4,9,1 1,13,14]  shows  that,  the  driving  fluid  mechanism  of  wing  rock  is  not 
unique  but  rather  is  greatly  dependent  on  the  model  geometry  and  configuration.  In  the  following, 
the  various  sources  of  wing  rock  will  be  discussed  based  on  the  experimental  and  computational 
results  provided  in  [4,7,22,29],  and  the  analysis  reported  in  [9,10,11,12,13,14,21,27,37]. 

1.2.1  Slender  Wing  Rock.  This  is  the  case  of  wing  rock  of  highly  swept-back  delta  wings 
and  limited  cases  of  blended  wing-bodies.  Slender  wing  rock  is  triggered  at  sufficiently  large 
angle  of  attack  by  an  asymmetric  leading  edge  (L.E.)  vortex  shedding.  Asymmetric  vortex 
shedding  could  be  the  result  of  flying  situations  such  as,  the  failure  of  the  pilot  to  recognize 
minute  incomplete  wing  leveling  when  he  pulls  up  on  the  stick.  Also,  flying  at  large  angle  of 
attack  in  asymmetric  flow  conditions  (typical  landing  in  a  cross-wind)  or  induced  lateral 
oscillations  due  to  unsteady  flow  over  the  wing  and/or  control  surfaces  could  lead  to  the  shedding 
of  an  asymmetric  vortex. 

To  completely  understand  the  fluid  mechanism  of  the  motion,  let  us  assume  a  slender 
delta  wing  (sweep  angle  >  76  degrees)  in  an  asymmetric  flow  condition  represented  by  a  slight 
positive  sideslip.  As  a  result,  the  L.E.  leeward  vortex  shifts  outboard,  and  the  L.E.  windward 
vortex  shifts  inboard  causing  the  wing  to  initially  roll  in  the  positive  roll  direction.  The  sudden 
roll  movement  causes  the  leeward  vortex  on  the  up-going  wing  to  compress  and  the  windward 
vortex  on  the  down-going  wing  to  stretch,  which  increase  the  initiated  rolling  moment.  As  the 
roll  angle  increases,  the  kinematic-coupling  between  the  angle  of  attack  and  the  sideslip  causes 
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the  effective  angle  of  attack  on  the  wing  to  decrease  and  the  effective  sideslip  angle  to  increase. 
The  increased  sideslip  on  the  wing  during  roll  causes  the  windward  vortex  on  the  down-going 
wing  to  move  inboard  and  toward  the  surface,  and  the  leeward  vortex  on  the  up-going  wing  to 
move  outboard  and  lifted  off.  The  convective  time  lag  associated  with  the  motion  of  the  lifted  off 
vortex  causes  the  right  wing  to  continue  dipping  until  the  lifted  vortex  takes  its  final  position. 
Then  the  lift  on  the  down-going  wing  is  increased,  while  that  on  the  up-going  wing  is  decreased. 
The  previous  effect  of  vortex  lift-off  together  with  the  reduction  of  vortex  strength  due  to  the 
decrease  of  the  effective  angle  of  attack  causes  the  wing  to  stop  at  a  finite  roll  angle  (the  limit 
cycle  amplitude),  and  then  reverses  its  motion.  As  the  wing  reverses  its  motion,  the  effective 
angle  of  sideslip  decreases  and  the  lifted  vortex  starts  to  reattach.  The  convective  time  lag  of  the 
vortex  motion  helps  the  rolling  moment  to  build  up  in  the  reversed  direction  until  the  vortex  is 
completely  reattached  to  the  leeward  side  of  the  left  wing.  What  remains  to  be  explained  is  how 
the  wing  overshoots  its  initial  equilibrium  (bank  zero).  This  is  due  to  the  hysteresis  associated 
with  the  structure  and  location  of  the  L.E.  vortices  over  the  two  wing  halves  during  the  oscillation 
as  shown  in  figure  1-1. 


Figure  1-1.  Vortex  structure  at  the  trailing  edge  during  rolling,  generated  by  the 

Vortex-Lattice  method  [29]. 


1-5 


1.2.2  Conventional-Wing  Rock.  This  is  the  wing  rock  associated  with  aircraft  of  straight 
or  moderately  swept  wings,  and  moderate  to  high  aspect  ratio.  In  other  words,  the  case  in  which 
the  aerodynamics  of  the  wing  is  dominated  by  the  two  dimensional  streamwise  flow,  and  wing 
rock  is  closely  related  to  dynamic  stall  [11]. 

The  fluid  mechanism  of  this  type  of  wing  rock  is  best  understood  by  considering  the 
analysis  of  plunging  oscillations  of  airfoils.  This  analysis  is  very  similar  to  that  of  the  well-known 
dynamic  stall  phenomenon  associated  with  a  nonzero  pitching  airfoil  [14].  Simply,  if  an  airfoil  is 
perturbed  in  plunge  when  flying  close  to  stall,  the  down-rolling  wing-half  will  experience  an 
upstream  movement  of  the  viscous  wall  effect,  which  promotes  separation  on  its  windward  side, 
and  hence  reduces  lift  over  the  static  one  (figure  1-2-A).  On  the  other  hand  the  up-rolling  wing 
half  will  experience  a  downstream  movement  of  the  wall  effect,  which  will  delay  the  separation 
on  its  leeward  side,  and  hence  increases  lift  over  the  static  one  (figure  1-2-B).  This  dynamic  stall 
effect  produces  the  driving  rolling  moment  for  the  wing.  As  the  rolling  angle  increases,  the 
kinematic-coupling  between  angle  of  attack  and  angle  of  sideslip  reduces  the  effective  angle  of 
attack.  At  a  sufficient  rolling  angle,  the  value  of  the  effective  angle  of  attack  goes  below  the  stall 
angle  of  attack,  which  causes  the  flow  to  reattach  on  the  windward  side  of  the  down-stroking 
airfoil.  This  generates  the  required  restoring  rolling  moment,  which  stops  the  roll  and  turns  it  in 
the  reverse  direction. 

At  high  subsonic  Mach  numbers,  the  plunging  effect  makes  little  difference  from  the 
previously  discussed  case.  The  flow  separation  in  this  case  is  driven  by  the  shock-boundary  layer 
interaction.  The  boundary  layer  transition  by  itself  causes  divergent  wing  bending  oscillations. 
The  moving  wall  effect  on  the  plunging  airfoil  of  the  bending  wing  will  promote  the  transition  on 
the  topside  during  up  stroke,  and  will  delay  the  transition  on  the  bottom  side  during  down  stroke, 
which  results  in  a  similar  effect  like  the  one  previously  discussed. 
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A- UPSTROKE 


Figure  1-2.  Wall  effect  on  an  airfoil  undergoing  plunging  motion  [14]. 

1.2.3  Wing-Body  Rock.  This  type  of  wing  rock  is  very  important  because,  it  is  obviously 
the  most  consistent  configuration  encountered  in  real  flying  vehicles.  Also,  experimental  results 
indicated  that  the  wing  rock  induced  by  fore-body  vortices  is  more  violent  than  the  one  induced 
by  the  leading  edge  of  a  flying  wing  alone  [4,13,21].  Typical  results  on  a  generic  fighter  model 
[4]  show  a  build-up  to  30-40  degrees  of  limit  cycle  amplitude  over  less  than  3  cycles.  This  is  very 
dangerous  in  the  sense  that  it  doesn’t  permit  enough  time  for  the  pilot  to  recover  before  the 
development  of  the  limit  cycle. 

According  to  the  flow  visualization  pictures  taken  during  the  wing  rock,  an  asymmetry¬ 
switching  mechanism  of  the  fore-body  vortices  was  encountered  in  the  same  fashion  to  that  of  the 
slender  delta  wing  leading  edge  vortices  [4,13]. 


The  following  discussion  on  the  effect  of  the  fore-body  cross-section  geometry  on  the 
wing  rock  characteristic  is  referred  to  Brandon  and  Nguyen  [4].  Different  fore-body  cross- 
sections  (vertical  ellipse,  hemispherical,  chine,  circular,  and  horizontal  ellipse)  are  mounted  on  a 
generic  fighter  model  (a  flat-plate  trapezoidal  wing  of  26  degrees  sweep,  flat-plate  horizontal  and 
vertical  tails,  and  fore-body  of  horizontal  ellipse  cross-section).  The  most  violent  oscillations 
were  found  for  the  horizontal  ellipse  (>30  degrees  amplitude),  the  moderate  oscillation  (20-25 
degrees  amplitude)  for  the  hemispherical  and  circular,  and  surprisingly  the  mild  oscillations  (<10 
degrees  amplitude)  for  the  vertical  ellipse  and  chine  (see  figure  1-3). 


Vertical  Hemispher-  Chine  Circular  Horizontal 

ellipse  leal  ellipse 


Figure  1-3.  Effect  of  various  fore-body  configurations  on  the  limit  cycle  amplitude  [4]. 

13  Overview 

The  remainder  of  this  study  consists  of  six  chapters.  In  Chapter-2,  the  previous  research  on 
wing  rock  simulations,  experimentally  and  computationally,  is  reviewed  and  presented.  In 
Chapter-3,  the  computational  model  and  the  methodology  adopted  in  the  solution  are  described  in 
detail.  Chapter-4  is  devoted  to  the  illustration  of  the  model  geometry  and  the  grid  generation 
process.  The  analysis  and  results  are  presented  in  Chapter-5,  with  the  conclusions  and 
recommendations  presented  separately  in  Chapter-6.  Also,  the  study  includes  four  appendices, 
appendix-A  and  B  contain  typical  input  and  output  for  the  PEGSUS  4.0  code  [51].  A  complete 
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derivation  for  the  Euler  boundary  conditions  is  found  in  appendix-C,  and  the  derivation  for  the 
computation  of  the  aerodynamic  forces  and  moments  is  presented  in  appendix-D. 
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Chapter  2-  Review  of  Previous  Work 


This  chapter  is  devoted  to  the  review  of  previous  work  on  the  experimental  and 
computational  simulation  of  wing  rock.  The  focus  in  the  review  is  centered  on  the  geometry  and 
configuration  of  the  tested  models,  the  utilized  degrees  of  freedom,  and  a  summary  of  the 
concluding  results. 

2.1  Experimental  Simulations 

The  earliest  observations  of  the  wing  rock  phenomena  were  reported  in  the  late  40’ s,  and 
since  then  numerous  studies  on  wing  rock  have  been  conducted  and  reported  in  the  literature. 
However,  the  understanding  of  the  fundamental  flow  mechanisms  that  drive  wing  rock  has 
remained  limited.  The  advancements  in  wind  tunnel  testing  techniques,  flow  visualization 
technology,  and  computational  fluid  dynamics  (CFD)  have  made  the  80’ s  and  90’ s  an  era  of 
recognized  intensive  experimental  and  computational  research  on  wing  rock. 

In  an  early  effort  to  clearly  understand  the  aerodynamic  mechanism  that  causes  slender  delta 
wings  to  rock,  the  NASA  Langley  Research  Center  (LaRC)  conducted  a  study  authored  by 
Nguyen,  Yip,  and  Chambers  in  1981  [42].  The  study  comprised  experiments  conducted  in  three 
low-speed  wind-tunnel  facilities  to  provide  four  kinds  of  testing,  standard  static-force,  forced- 
oscillation,  rotary,  and  free-to-roll  tests.  An  80°-sweep  flat-plate  delta  wing,  constructed  of  wood 
and  sharply  beveled  at  the  leading  edges  was  tested.  An  important  feature  of  the  model  was  that 
the  roll-axis  laid  in  the  plane  of  symmetry  parallel  to  the  x-body  axis  but  shifted  below  by  2 
inches.  The  research  provided  groundbreaking  results  for  the  slender  delta-wing  rock  with  single 
degree-of-freedom  (DoF)  in  roll.  The  model  started  to  rock  at  angle  of  attack  (AOA)  of  25 
degrees,  while  vortex  bursting  was  not  observed  below  30  degrees.  These  early  observations 
demonstrated  that  slender  wing  rock  is  usually  triggered  by  leading-edge  vortex  asymmetry, 
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while  vortex  burst  (if  encountered)  affects  the  value  of  the  amplitude  and  frequency  of  the  limit 
cycle.  Static  tests  for  various  sideslip  angles  showed  that  the  vortex  burst  location  on  the 
windward  side  moves  toward  the  leading  edge  and  the  leeward  vortex  core  is  lifted-off  the  wing. 
That  is,  sideslip  is  responsible  for  driving  one  of  the  two  major  aerodynamic  mechanisms  that 
control  the  wing  rock  motion.  Dynamically  speaking,  one  can  say  that  wing  rock  is  sustained  by 
the  dependence  of  the  aerodynamic  damping  in  roll  on  sideslip  such  that  it  is  positive  (unstable) 
for  small  sideslip  during  the  onset  of  wing  rock  and  negative  (stable)  for  large  sideslip  during  the 
fully-developed  wing  rock.  This  idea  was  actually  the  initial  motivation  for  the  present  research 
to  determine  the  significance  of  including  the  sideslip  DoF  on  the  limit  cycle  characteristics  of 
wing  rock. 

Subsequently  in  1982,  Levin  and  Katz  [31]  conducted  a  similar  experiment.  The  model 
geometry  was  very  much  similar  to  Nguyen’s  model  except  for  the  roll  axis,  which  coincided 
with  the  body  axis.  The  shifting  of  the  roll  axis  in  the  Nguyen’s  experiment  was  a  source  of 
slightly  increased  sideslip  (one  can  easily  prove  that,  A(3  =  pd/U where  p  is  the  roll  rate,  d  is  the 
shifted  distance,  and  U„  is  the  free-stream  velocity).  The  authors  were  able  to  construct 
boundaries  for  the  slender  wing  rock,  which  were  correlated  to  the  initial  combination  of  AOA 
and  sideslip  [31:36]. 

In  1986  NASA  LaRC  sponsored  research  aimed  at  investigating  the  effect  of  forebody 
geometry  on  the  static  and  dynamic  stability  of  aircraft  at  high  angles  of  attack.  The  research  was 
conducted  by  Brandon  and  Nguyen  [4],  and  the  results  were  previously  demonstrated  in 
Subsection  1.2.3. 

In  1986  Hall  and  Del  Frate  [21]  conducted  a  water-tunnel  study  on  the  interaction  between 
the  forebody  and  wing  primary  vortices.  The  authors  selected  three  common  fuselage  cross- 
sections,  circular,  chine  of  90°-included  angle,  and  chine  with  7.5°-included  angle  that  produce 
weak,  medium,  and  strong  forebody  vortex  strengths  respectively.  This  analysis  of  the  effect  of 
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three  forebody  cross-sections  concluded  that  the  stronger  the  forebody  vortex,  the  more  energy  is 
imparted  to  the  wing  vortex  system,  and  consequently  the  more  resistance  to  vortex  burst.  The 
results  for  the  three  forebodies  mounted  on  a  55°-sweep  cropped  delta  wing  model  in  zero 
sideslip  showed  a  burst  onset  of  22.5°  in  AOA  for  the  circular,  40°  for  the  90°-chine,  and  45°  for 
the  7.5°-chine.  In  sideslip,  large  differences  in  the  burst  locations  between  the  windward  and 
leeward  wing  sides  were  clearly  distinguished,  which  in  turn  increases  the  asymmetry  tendency 
in  sideslip. 

Jun  and  Nelson  [24]  in  1987  conducted  experiments  on  an  80°-sweep  delta  wing  similar  to 
the  experiments  found  in  [31,42],  With  only  slight  differences  in  the  model  installation,  which 
was  mounted  on  a  sting  that  permits  free  roll  of  the  model  about  a  bearing  system  located  in  the 
sting  support.  The  experimental  results  confirmed  the  early  results  in  [31,42]  that  the  dynamic 
vortex  trajectories  differ  from  the  static  trajectories  and  depend  on  the  direction  of  roll. 
Differences  were  also  observed  in  the  location  of  vortex  breakdown  between  the  static  and 
dynamic  cases. 

In  1993,  Straka  and  Hemsch  [50]  conducted  an  experiment  in  the  NASA  LaRC  water-tunnel 
to  determine  the  effect  of  a  circular-section  forebody  on  the  location  and  onset  of  vortex  burst  on 
a  69.3°-sweep  delta  wing.  The  variation  in  vortex  burst  location  with  angle  of  attack  at  zero 
sideslip  was  recorded  for  the  wing  alone,  and  the  wing-fuselage  combination.  It  was  concluded 
that,  adding  the  body  to  the  delta  wing  adversely  affected  the  onset  and  the  chordwise 
progression  of  vortex  burst  on  the  wing. 

Preliminary  prediction  of  aircraft  stall,  post-stall,  and  spin  characteristics  based  on  classical 
prediction  methods  or  systematic  wind  tunnel  testing  are  limited  due  to  the  complexity  of  the 
aerodynamics.  The  most  reliable  source  of  information  about  those  critical  maneuvers  prior  to 
actual  flight  tests  is  the  unconstrained  testing  of  a  properly  dynamically-scaled  model.  The 
principal  techniques  were  the  wind-tunnel  free-flight  technique  for  stall  and  departure,  the  spin- 
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tunnel  technique  for  spin  and  recovery,  and  the  radio-controlled  drop-model  technique  (RCDMT) 
for  stall  and  post-stall  motions.  In  1987  Fratello  et  al  [15]  conducted  research  at  NASA  LaRC  to 
study  the  high-AOA  characteristics  for  a  20%  dynamically-scaled  X-29A  model,  utilizing  the 
unconstrained  RCDMT.  The  wind-tunnel  tests  had  predicted  a  25-degrees  angle  of  attack  for  the 
onset  of  wing  rock.  The  results  obtained  from  the  RCDMT  verified  an  earlier  onset  at  20°,  and 
with  further  increase  of  angle  of  attack  the  wing  rock  motion  grew  rapidly,  with  the  amplitude 
doubled  in  only  one  oscillation.  At  an  average  angle  of  attack  of  32°,  the  motion  diverged  into  a 
violent  snap  roll  departure  of  peak  roll  rate  over  (130  deg/s). 

In  summary,  the  previous  experimental  simulations  on  wing  rock  have  been  limited  to  a 
single  DoF  in  roll.  The  bearing  mounting  and  friction  have  affected  the  predicted  limit  cycle 
characteristics  as  was  previously  reported  by  Levin  and  Katz  [31].  Only  the  RCDMT  can  provide 
a  realistic  simulation;  however  it  is  still  considered  expensive  and  technically  difficult,  and  hence 
of  limited  use. 

2.2  Computational  Simulations 

The  rapid  developments  in  computer  speed,  architecture,  and  memory  in  the  early  80’s  had  a 
great  impact  on  the  improvement  of  aerodynamic  numerical  solutions.  The  improvements  in  the 
aerodynamic  numerical  solutions  have  been  complemented  by  parallel  progress  in  grid 
generation  and  flow  visualization  techniques.  The  whole  package,  the  grid  generation,  flow 
solver,  and  flowfield  visualization  utilities,  have  become  the  computational  tunnel  for  present- 
day  aerodynamic  research.  Most  of  the  developed  computational  codes  have  been  successfully 
validated  by  comparison  with  numerous  experimental  data. 

Computational  simulation  of  wing  rock  has  been  encouraged  due  to  the  difficulty  of  model 
installation  for  free  motion  in  the  tunnel,  and  for  better  flowfield  visualization.  In  the  following 
review,  the  computational  work  on  wing  rock  simulation  will  be  reviewed  and  presented. 
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In  1982  Levin  D.  and  Katz  J.  [32]  introduced  a  modified  version  of  the  classical  Vortex 
Lattice  method  that  is  capable  of  solving  both,  the  steady  and  unsteady  aerodynamics  as  well, 
with  an  application  for  the  simulation  of  leading  edge  vortices  over  delta  wings  at  high  AOA. 

Shortly  afterwards  Konstadinopoulos,  Mook,  and  Nayfeh  [29]  coupled  the  unsteady  vortex- 
lattice  method  with  the  rigid  body  dynamics  in  roll  to  numerically  simulate  the  previous 
experiments  for  Nguyen  [42]  and  Levin  [31].  The  results  showed  good  agreement  with  the 
corresponding  experimental  data. 

In  1989  Katz  J.  [26]  introduced  a  new  application  for  the  Unsteady  Vortex  Lattice  method 
by  computing  the  conical  rotary  loads  on  a  standard  dynamic  model.  The  author  was  able  to 
demonstrate  excellent  agreement  with  the  experimental  data,  at  low  Mach  number  (M  =  .15)  and 
moderate  angles  of  attack  (AOA  <15  degrees),  with  extremely  low  computational  cost. 

The  limited  capabilities  of  the  early  supercomputers  nudged  the  first  CFD  simulations  of 
wing  rock  to  less-costly  computational  solutions  such  as  the  Conical  Euler  and  Euler  schemes. 
One  of  the  earliest  works  utilizing  these  schemes  was  due  to  Lee  and  Batina  [30]  in  1991  from 
NASA  LaRC.  The  simulation  has  been  done  for  a  slender  delta-wing  configuration  (leading  edge 
sweep  of  75°)  at  Mach  number  1.2.  The  supersonic  freestream  Mach  number  allowed  the 
implementation  of  the  conical  flow  similarity,  and  hence  reduced  the  problem  from  three 
dimensions  to  two  dimensions,  which  considerably  reduced  the  computational  cost.  The  authors 
utilized  the  cell-centered  finite-volume  algorithm  for  solving  on  an  unstructured  grid  made  up  of 
triangles.  An  explicit  artificial  dissipation  of  second  and  fourth  order  was  added  to  prevent 
oscillations  near  shock  waves  and  damp  out  high-frequency  uncoupled  error  modes.  The  rigid- 
body  roll  equation  has  been  explicitly  coupled  with  the  flow  governing  equations.  The  time 
derivatives  in  the  rigid  body  equation  were  approximated  by  a  2nd-order  accurate  backward  finite- 
difference.  The  number  of  elements  used  in  this  study  was  8299,  which  had  a  total  of  4226 
nodes,  with  refinement  on  the  leeward  side  of  the  wing  where  the  dominant  flow  features  were 
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expected  to  occur.  The  wing  rock  was  excited  by  an  initial  roll  rate  at  angle  of  attack  of  30°,  and 
a  limit  cycle  of  amplitude  38°  was  observed,  which  agreed  well  with  the  earlier  experiments. 

In  1992,  Kandil  and  Salman  [25]  took  a  step  forward  by  coupling  the  unsteady  Euler 
equations  with  the  rigid  body  dynamics  in  roll.  They  simulated  the  wing  rock  oscillations  and 
divergence  for  an  80°-sweep  delta  wing  at  30°  AOA  and  Mach  number  of  0.3.  The  unsteady 
Euler  equations  was  used  to  simulate  the  fluid  dynamics,  while  the  four-stage  Runge-Kutta 
scheme  was  used  to  solve  the  rigid-body  dynamic  equations.  The  flow  governing  equations  were 
solved  in  the  moving  frame  of  reference,  which  adds  an  extra  computational  cost  for  the  source 
terms.  However,  the  use  of  this  approach  saves  regenerating  the  grid  and  recalculating  the  values 
of  the  transformation  metrics  after  each  time  step.  The  implicit,  approximately-factored, 
centrally-differenced,  finite-volume  scheme  was  employed,  with  2nd  and  4th-order  explicit 
dissipation  terms.  Implicit  2nd-order  dissipation  terms  were  also  added  to  enhance  the  scheme 
stability.  The  grid  consisted  of  32x32x48  grid  points  in  the  axial,  normal,  and  wrap-around 
directions. 

In  1993  Chaderjian  and  Schiff  [6]  presented  an  early  application  of  a  generalized  full 
Navier-Stokes  based-code.  They  conducted  static  tests  to  investigate  the  nonlinear  variation  of 
the  rolling  moment  coefficient  with  angle  of  attack  for  a  65°-sweep  delta  wing  at  30°  AOA, 
3.67xl06  Reynolds  number,  and  a  Mach  number  of  0.27.  The  three-dimensional,  Reynolds- 
averaged,  Navier-Stokes  (RANS)  equations  were  used  to  model  the  viscous  flow.  The  relatively 
high  Reynolds  number  used  in  this  study  allowed  the  application  of  thin-layer  approximation. 
The  perfect  gas  law,  Sutherland’s  viscosity  law,  and  a  turbulence  model  were  used  to  complete 
the  RANS  system  of  equations.  The  governing  equations  were  transformed  to  body-fitted 
curvilinear  coordinates  using  general  coordinate  transformation.  The  numerical  procedure 
comprised  of  the  implicit  Beam  and  Warming  algorithm  [3],  approximately-factored,  and 
centrally-differenced.  A  three-dimensional  hyperbolic  grid  generator  was  used  to  generate 
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67x49x208  nodes  in  the  axial,  normal,  and  circumferential  directions,  a  total  of  about  700,000 


nodes.  The  run  of  this  code  for  one  static  roll  angle  was  reported  to  require  48  hours  on  the  C-90 
computer  to  reach  a  steady  state. 

In  an  effort  to  validate  the  code  developed  in  [6]  for  dynamic  testing,  Chaderjian  [5]  in  1994 
simulated  the  forced  large-amplitude,  high  frequency  roll  oscillations.  The  RANS  equations  were 
written  in  the  inertial  frame  of  reference  and  transformed  to  the  body-fitted  curvilinear  frame  of 
reference  via  the  generalized  time-dependant  coordinate  transformation.  The  comparison  of  the 
results  with  the  corresponding  experimental  data  showed  a  very  good  quantitative  agreement. 
The  computational  cost  of  the  developed  code  to  advance  one  time  step  on  the  Cray  Y-MP  is  18 
jus/grid-point/time-step.  One  cycle  of  oscillations  took  15,000  time  steps,  which  requires  50  hours 
of  CPU-time  for  a  total  of  700,000  grid  points.  Later  in  1996,  Chaderjian  and  Schiff  [7]  coupled 
the  rigid-body  dynamic  equations  in  roll  with  the  RANS  equations  and  simulated  the  free  roll 
motion  of  wing  rock  for  the  same  previous  model.  The  resulted  computational  free  roll  motions 
were  reported  to  damp  out  more  rapidly  compared  with  the  experimental  simulations,  but  the 
computed  and  experimented  frequencies  were  in  good  agreement. 

In  1993,  Gordnier  and  Visbal  [19]  presented  an  application  for  th tfdl3cli  code,  which  has 
gone  through  a  similar  validation  procedure  from  static  to  forced  roll  tests.  The  code  is  very 
similar  to  the  one  developed  by  Chaderjian  and  Schiff  [6],  except  for  a  few  differences.  Unlike 
the  previous  code,  which  utilizes  the  thin-layer  version  of  the  Navier-Stokes  equations,  th tfdl3di 
utilizes  the  unsteady,  three-dimensional,  full  mass-averaged  Navier-Stokes  equations,  which 
extends  its  application  to  low  Reynolds  numbers  as  well.  An  additional  feature  in  th efdl3di  code 
is  the  implementation  of  the  sub-iteration  algorithm  to  reduce  the  various  errors  due  to  the 
factorization,  the  lag  in  applying  the  boundary  conditions,  and  the  inherent  lag  in  the  solution  of 
the  coupled  equations.  The  sub-iteration  algorithm  has  greatly  improved  the  stability  and 
accuracy  of  the  computational  scheme.  The  code  has  been  validated  by  experimental  data  of  two 
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previous  experiments,  for  the  80°  and  65°  delta  wing  models.  Both,  static  and  forced-roll 
simulations,  were  conducted  to  demonstrate  the  variations  in  the  vortex  structure  and  strength 
due  to  the  dynamic  interaction  during  a  constant  roll  rate  motion.  A  useful  discussion  on  the 
vortex  dynamics  has  been  introduced  to  explain  the  reduction  in  the  restoring  moment  in  the 
dynamic  case  (the  presented  analysis  on  slender  wing  rock  in  subsection  1.2.1  is  mostly  referred 
to  this  paper).  Gordnier  [16]  has  analyzed  the  effect  of  roll  rate  and  initial  acceleration  on  the 
vortex  dynamics  during  a  forced-roll  maneuver  for  the  80°-sweep  delta  wing  model.  In  1997, 
Gordnier  [18]  presented  static  roll  tests  for  the  same  model  used  by  Chaderjian  [5,6,7]  but  at 
lower  Reynolds  and  Mach  numbers  to  demonstrate  the  ability  of  the  code  to  simulate  complex 
viscous  flows.  Also,  Gordnier  [17]  has  demonstrated  the  ability  of  the  fdl3di-code  to  simulate  the 
spiral  vortex  breakdown  and  the  transition  from  spiral  to  bubble-type  breakdown  with  very  good 
qualitative  agreement  with  the  experimental  data. 

In  an  effort  to  model  the  nonlinear,  unsteady  aerodynamics  for  the  rolling  65°  delta  wing, 
Tromp  [52]  has  simulated  the  constant-rate  conical  motion,  constant-roll  rate,  and  the  barrel-roll 
maneuver.  The  numerical  simulation  has  been  done  using  thefdl3di-code  for  zero  and  30  degrees 
AOA.  The  numerically  computed  loads  have  shown  an  excellent  agreement  with  the 
experimental  data. 

In  summary,  it  has  been  shown  in  the  previous  review  that  the  computational  wing  rock 
simulations  were  commonly  conducted  for  a  single  DoF  in  roll  and  applied  to  either  flying  delta 
wing  or  blended  wing-body  (without  fore-body)  configurations. 
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Chapter  3  -  Computational  Model 


3.1  Flow  Aerodynamics 

All  the  computational  experiments  in  this  study  are  conducted  at  a  Mach  number  of  0.4  and 
pressure  altitude  of  2000  meters.  The  corresponding  Reynolds  number  for  this  flight  regime, 
based  on  the  root  chord  length,  is  significantly  high  (ReM  =  2.472xl09),  which  motivated  the  use 
of  the  Euler  equations  to  reduce  the  computational  cost.  The  Euler  equations  have  been  used 
widely  to  simulate  vortical  flow  aerodynamics  on  highly  swept-back  or  delta-wing  configurations 
at  moderate-to-high  values  of  Reynolds  number  [1,25,30,38,47].  The  Advisory  Group  for 
Aerospace  Research  &  Development  (AGARD)  [1:162]  in  its  review  on  computational 
aerodynamics  based  on  the  Euler  equations  has  reported  that  the  Euler  solutions  have  the  ability 
to  adequately  simulate  vortical  flow  aerodynamics  and  leading  edge  vortex  breakdown  on  delta¬ 
wing  configurations.  Also,  Yu  et  al  [55]  have  reported  excellent  agreement  with  experimental 
data  for  the  Euler  solution  of  a  complete  commercial  transporter  in  transonic  flow  conditions  and 
sideslip.  A  modified  version  of  the  fdl3di  code,  provided  by  the  Air  Force  Research  Laboratory 
(. AFRL ),  Wright-Patterson  Air  Force  Base,  Ohio,  has  been  utilized  to  provide  the  finite  difference 
modeling  and  solution  of  the  Euler  flow  governing  equations. 

3.1.1  Flow  Governing  Equations.  The  unsteady  three-dimensional  compressible  Euler 
equations  in  strong  conservation  form  have  been  used.  The  equations  have  been  written  in  a  fixed 
inertial  frame  of  reference  and  transformed  to  the  computational  domain  using  a  generalized 
time-dependant  transformation  (£,77  ,t).  The  flow  variables  are  nondimensionalized  as  follows: 

U*  —  U/U„ 
v*  =  v/£/„ 
w*  =  w/U „ 

P*  =  p/(  p„  uJ) 
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p  =  p/p~ 
T*  =  7771 


e*=  e/tA,2 

All  coordinate  lengths  have  been  normalized  by  the  root  chord-length  of  the  delta  wing,  Cr. 
x*  =  x/Cr 
y*  =  y/Cr 

Z*  =  Z/Cr 


Where  x  =  y  =  z  =  0  is  at  the  nose. 

The  final  form  of  the  transformed  nondimensionalized  equations  in  strong  conservation  form  can 
be  written  as  follows: 


BQ  dF_  dG  dH_ 

dt  +  dE,  drj  dC 


(3-1) 


Where: 


The  vector  of  dependent  variables  is  given  by: 
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*  * 

p  u 

1 

* 

_ I 

*  * 

p  w 

p*  uU  +  ZXP* 

G  =  — 

J 

*  *T  T  «  ^ 

p  u  V  +rixP 

H- 1 

J 

p*u*W+ZxP* 

p*v*u  +£,yP* 

*  *77  n* 

P  V  V  +T]yP 

PVw+CyP * 

pw'u+^p* 

J 

p*  w*V  +t]zP* 

J 

p*w*W+CzP* 

p*  e*  U  +  P*U' 

p*e*V+P*V\ 

p*e*W  +  P*W\ 

u=Zl+Zxu+Z,v-+SS 

=  l+u' 


3-2 
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=  V,+V' 
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=  C.+W' 
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The  Jacobian  of  the  transformation,  J,  is  given  by: 


J  = 


777“*  *  +1  *7  *  *  *  H  *7  *  *  *  *  \  ,  *"7  *  *  * 

3(*  ,y  ,z  )  x(\ynz(-y(zn)-Xi\yfzc-ycZf)+x(\y(zn-y,lZt) 

The  perfect  gas  equation  closes  the  system  of  equations  (3-1)  as  follows: 

PT * 


P  =■ 


yM' 


(3-2) 


3.2.2  Boundary  Conditions.  The  following  boundary  conditions  are  applied  to  the  flow 
governing  equations: 

Characteristic  boundary  conditions  [54]  are  applied  at  the  far-field  and  upstream  boundaries. 
The  values  of  the  flow  variables  at  the  exit  boundary  are  extrapolated  from  the  interior  by 
first-order  extrapolation. 

Periodic  boundary  conditions  have  been  applied  for  the  points  of  symmetry  in  the  wrap¬ 
around  direction. 

On  the  body  surface  the  classical  Euler  boundary  conditions  have  been  applied: 

U„=Ub-n 


dP_ 

dn 


-pab  -n 
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Where  Un  is  the  flow  velocity  component  in  the  direction  normal  to  the  body  surface.  U b  and  ab 
are  the  velocity  and  acceleration  at  the  surface  grid  point  location  rb ,  given  by: 
Ub=UC!i+6jx(rb-rJ 

+^Xfe  -rJ-H0x[<S>x(i-  -rj] 

3.1.3  Numerical  Solution  of  the  Flow  Governing  Equations.  The  governing 
equations  are  solved  numerically  using  the  implicit,  approximately-factorized  finite  difference 
scheme  of  Beam  and  Warming  [3]  in  the  diagonalized  form  of  Pulliam  and  Chaussee  [46].  The 
equations  are  differenced  using  second-order  accurate  differencing  for  temporal  and  spatial 
derivatives.  A  blend  of  nonlinear  second  and  fourth  order  spectral  damping  derivatives  is 
employed.  The  boundary  conditions  are  explicitly  applied,  and  a  sub-iteration  algorithm  is 
employed  to  reduce  the  effect  of  the  inherent  time  lag  in  applying  the  boundary  conditions  and 
reduce  the  factorization  error. 

The  numerical  implementation  of  the  boundary  conditions  on  the  body  surface  is  derived  and 
explained  in  detail  in  Appendix-C.  Also,  the  numerical  computation  of  the  aerodynamics  forces 
and  moments  on  the  body  surface  is  presented  in  Appendix-D. 

3.2  Rigid  Body  Dynamics 

The  relevant  degrees  of  freedom  in  this  study  are  the  roll,  sideslip,  and  vertical  motion.  The 
rolling  equation  is  written  in  the  body-axes  frame  of  reference  to  keep  the  roll-axis  moment  of 
inertia  constant  throughout  the  entire  motion,  while  the  sideslip  and  vertical  motion  equations  are 
written  in  the  inertial  frame  of  reference. 

3.2.1  Rigid  Body  Dynamic  Equations.  The  rolling  equation  is  a  second-order 
autonomous  ordinary  differential  equation  in  time. 
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Where  Lh  is  the  body-axis  rolling  moment,  and  Ix  is  the  moment  of  inertia  about  the  body-fixed  x- 
axis. 

The  sideslip  and  vertical  equations  are  written  in  the  inertial  frame  of  reference  since  there  is  no 
need  to  be  solved  in  the  body-axes  frame  of  reference.  The  equations  are  second-order 
autonomous  ordinary  differential  equations  in  time. 
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Where  Fy  and  Fz  are  the  side  and  vertical  force  components  computed  in  the  y  and  z  inertial 
directions  respectively. 

A  first  order  state-space  representation  of  (3-3),  (3-4),  and  (3-5)  is  given  by: 
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Introduce  the  non-dimensional  forms  of  the  roll  rate,  time,  rolling  moment  of  inertia,  mass,  and 
gravity  acceleration  respectively  as  follows: 
pb 
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The  first-order  system  (equation  3-6)  is  given  in  non-dimensional  form  as  follows: 
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(3-7-b) 

(3-7-c) 

(3-7-d) 

(3-7-e) 


(3-7-f) 
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3.2.2  Numerical  Solution  of  the  Rigid  Body  Dynamic  Equations.  The  non- 
dimensional  first-order  system  (3-7)  has  been  differenced  in  time  using  second-order  accurate 
backward  difference.  Thus,  the  finite  difference  solution  at  time  level  (n)  is  given  as  follows: 
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(3-8-a) 

(3-8-b) 

(3-8-c) 

(3-8-d) 

(3-8-e) 


(3-8-f) 


Where  the  superscripts  between  brackets  denote  the  corresponding  time  level. 


3.3  Coupling  the  Flow  Governing  Equations  and  the  Rigid  Body  Dynamic  Equations 
The  rigid  body  equations  are  coupled  explicitly  with  the  flow  governing  equations.  As  seen 
in  equations  (3-8-a,  c,  and  e),  the  aerodynamic  forces  and  moments  are  computed  at  a  time  level 
(n-1).  This  inherent  lag  in  the  solution  of  the  coupled  equations  has  been  addressed  through  a  sub¬ 
iteration  algorithm.  The  original  fdlSdi  code  already  employed  a  Newton-like  sub-iteration 
algorithm  on  the  flow  governing  equations,  which  was  extended  to  operate  on  the  coupled  set  of 
equations.  Figure  3-1  depicts  a  simple  flowchart  demonstrating  the  application  of  the  sub¬ 
iteration  approach  on  the  coupled  set  of  equations. 
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3.4  Implementation  of  the  Body  Motion 

The  flowfield  interactions  due  to  the  body  dynamics  are  implemented  in  the  modified  fdl3 di¬ 


code  as  follows: 


After  solving  the  rigid-body  equations  the  new  position  of  the  individual  grids,  in  the  inertial 


coordinate  system,  is  computed  by  the  following  transformation: 
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Where,  the  superscript  (0)  denotes  the  initial  grid  position  at  time  level  (n  =  0),  which 
corresponds  to  the  initial  states  (<()(0)  =  0,  ycg*(0)  =  0,  and  zcg  (0)  =  0). 

Compute  the  new  values  of  the  space  and  time  metrics  of  the  generalized  coordinate 
transformation. 

Compute  the  velocities  and  accelerations  on  the  body  surface  grid  points,  which  are  needed 
for  the  application  of  boundary  conditions. 


3.5  Computer  Resources 

All  the  computational  simulations  presented  in  this  study  has  been  run  on  the  CRAY  T-94 
supercomputer  at  the  Ohio  Supercomputer  Center  (OSC),  Columbus,  Ohio,  as  part  of  the  granted 
CPU-Hours  for  project  number  PIS  198. 

The  CRAY  T-94  computer  system  at  the  OSC  is  a  powerful  general-purpose  supercomputer 
that  features  4  high-speed  (450  MHz)  processors,  each  with  a  peak  performance  of  approximately 
2  billion  floating  point  operations  per  second  (2  GFLOPS). 
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Chapter  4  -  Grid  Generation 


4.1  The  Model 

The  tested  model  in  this  study  is  a  5%  dynamically-scaled  model  for  a  generic  fighter 
consisting  of  fore-body,  blended  wing-body,  and  vertical  tail  (see  figures  4-1  and  4-2).  The  fore¬ 
body  cross  section  is  a  chine  of  90°-included  angle  with  a  sharp  nose  of  15°  apex-half  angle.  The 
wing  is  a  thin  flat  plate  cropped  delta  wing  of  65°  sweep  with  the  leading  and  trailing  edges 
beveled  at  45°.  The  vertical  tail  has  a  trapezoidal  plane-form  with  beveled  edges,  and  a  finite 
thickness.  The  geometry  of  the  body,  wing,  and  tail  is  chosen  for  simplicity,  but  in  fact  is  adopted 
from  realistic  present-day  fighters  such  as  the  Dassault-Breguet  Mirage  [23]. 


*  All  dimensions  are  in  cm,  unless  otherwise  stated. 


Figure  4-1.  Three-view  drawing  of  the  generic  fighter  model. 
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Figure  4-2.  Three-dimensional  view  of  the  model. 

Originally,  the  model  mass  was  to  be  estimated  from  the  mass  of  the  Mirage  2000  [23], 
(7440.3  Kg,  empty).  Due  to  the  unavailability  of  the  values  of  the  moment  of  inertia  for  the 
Mirage,  the  mass  and  longitudinal  moment  of  inertia  for  the  model  were  estimated  instead  from 
that  of  the  A-4  fighter,  which  has  a  similar  size  and  shape  to  the  Mirage.  The  inertia  properties  of 
the  A-4  fighter  [40:405]  are  given  by: 
m  =  7910  [Kg],  and  Ix  =  10968  [kg-m2]. 

So,  the  non-dimensional  inertia  parameters  of  the  model  (refer  to  subsection  3.2.1)  are  given  as 
follows: 

m*=  35.3,  and  I*  =  0.684 

4.2  Meshing  the  Model 

The  accuracy  of  a  CFD-simulation  is  in  general  governed  by  both,  the  quality  of  the  mesh 
that  resolves  the  relevant  model  geometry  and  physics,  and  the  time/spatial  accuracy  of  the 
computational  scheme.  This  emphasizes  the  role  of  the  grid  generation  process  in  order  to  get 
high  quality  simulations.  Before  starting  the  grid  generation  process,  four  key  elements  should  be 
carefully  selected: 
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•  The  grid  structure  (structured  vs.  unstructured). 

•  The  grid  generation  technique  (Single  vs.  Multi-overset  mesh  techniques). 

•  The  grid  topology  for  each  individual  mesh  (H-O,  H-H,  or  C-O. . .). 

•  The  domain  dimensions. 

In  the  remaining  part  of  this  section,  the  author’s  choice  for  these  elements  is  discussed  and 
justified. 

The  utilized  code  in  this  study  to  solve  the  flow  governing  equations  is  the  fdl3di-c ode, 
which  is  written  for  structured  grids,  leaving  no  choice  for  grid  structure. 

In  order  to  justify  the  selection  for  the  meshing  technique,  one  needs  to  differentiate  among 
the  available  alternatives.  Three  alternatives  are  possible:  either  a  single  O  or  H-grid,  or  multi¬ 
overset  grids.  Two  factors  were  taken  into  consideration  in  this  investigation,  the  total  number  of 
grid  points,  and  the  quality  of  the  generated  mesh.  In  the  case  of  a  single  O-grid,  the  number  of 
grid  points  must  be  increased  in  the  wrap-around  direction  at  the  location  of  the  vertical  fin  to 
resolve  its  geometry.  This  increase  of  grid  points  must  be  made  all  the  way  along  the  axial 
direction,  resulting  in  unnecessary  augmentation  of  the  total  number  of  grid  points.  Moreover,  the 
generation  of  a  single  O-structured  grid  for  the  whole  model  including  the  tail  was  found  to  be 
difficult  if  not  impossible.  For  the  case  of  a  single  H-grid,  the  clustering  of  grid  lines  to  resolve 
the  model  geometry  will  propagate  all  over  the  computational  domain,  resulting  in  a  significant 
increase  in  the  total  number  of  grid  points  (possibly  greater  than  that  for  the  case  of  a  single  O- 
grid).  The  quality  of  the  resulting  mesh  in  the  case  of  a  single  H-grid  may  exceed  the  one 
obtained  by  the  O-grid.  However,  the  cut  in  the  H-grid  made  by  the  solid  body  surface  requires 
modifying  the  computational  stencils  at  these  points  to  be  one-sided  everywhere  in  the  code.  On 
the  other  hand,  the  utilization  of  the  multi-overset  mesh  technique  allows  independent  meshing  of 
the  sub-domains  of  interest,  which  are  overset  and  mutually  communicate  across  the  boundaries. 
In  addition  to  the  flexibility  obtained  by  the  utilization  of  the  multi-overset  technique,  the  overall 
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number  of  grid  points  is  highly  reduced  compared  with  that  for  the  case  of  a  single  mesh.  For  all 
the  previous  reasons,  the  multi-overset  mesh  technique  was  found  to  be  best  suited  for  our 
particular  model.  The  model  is  subdivided  into  two  meshes,  one  for  the  wing-body  and  the  other 
for  the  vertical  fin.  The  fin  mesh  is  embedded  in  the  wing-body  mesh,  and  a  final  composite  mesh 
is  created. 

4.2.1  The  Wing-Body  Mesh.  The  axi-symmetric  geometry  of  the  wing-body  and  the 
conical  shape  of  the  fore-body  motivate  the  utilization  of  a  C-0  topology  for  this  mesh.  Also,  it 
becomes  relatively  easier  to  apply  the  boundary  conditions  on  the  body  surface  in  case  of  a  C-0 
grid.  The  C-0  topology  has  been  extensively  used  in  the  literature  to  mesh  blended  wing-body 
configurations  [5,6,7,18,25,52].  A  structured  mesh  has  been  generated  using  Gridgen  ver.  13 
[44];  an  initial  grid  is  generated  for  a  quarter  of  the  domain  using  polar  transfinite  interpolation. 

An  elliptic  partial  differential  equation  ( EPDE)  solver  was  run  to  improve  the  quality  of  the 
grid,  and  control  the  grid  lines  near/on  the  body  surface.  A  Thomas-Middlecoff  control  function 
[44]  was  applied  to  map  the  clustering  of  the  grid  points  on  the  boundary  to  the  interior  domain. 
A  Steger-Sorenson  control  function  [44]  was  applied  on  the  body  surface  to  keep  the 
orthogonality,  and  hence  simplifies  the  application  of  boundary  conditions.  Due  to  the  axi- 
symmetric  geometry  of  the  wing-body  section,  only  a  quarter  of  the  domain  is  generated  by 
Gridgen,  and  the  other  three  quarters  are  generated  by  reflection  using  a  FORTRAN  routine.  The 
advantage  of  doing  this  is  to  enforce  an  exact  symmetry  on  the  distribution  of  mesh  points,  which 
may  be  slightly  violated  if  the  PDE  solver  has  been  run  for  the  whole  domain.  Three  grids  of 
different  dimensions  have  been  generated  (Table  4-1). 

The  domain  consists  of  a  circular  cylinder  connected  to  a  hemisphere  (figure  4-4).  The 
length  of  the  circular  cylinder  is  three  characteristic  lengths  (wing  root  chord  length  Cr),  and  has  a 
radius  of  1.5  characteristic  lengths  (figure  44).  The  domain  shape  and  dimensions  have  been 
previously  utilized  in  most  meshes  of  C-0  topology  [18,25]. 
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Figure  4-3.  Section  in  the  wing-body  mesh  (coarse)  demonstrates  the  topology. 


Table  4-1.  Dimensions  of  the  grids. 


Grid  Type 

Dimensions 

Total 

No.  of 
Points 

(axial) 

r|  (normal) 

£  (wrap-around) 

Coarse 

79 

35 

65 

179,725 

Medium 

117 

53 

97 

601,497 

Fine 

120 

70 

105 

882,000 
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4.2.2  The  Vertical  Fin  Mesh .  The  vertical  fin  mesh  consists  of  the  fin  and  part  of  the 
fuselage  in  order  to  resolve  the  area  at  the  fin-fuselage  connection  (see  figure  4-5).  An  H-H  grid 
is  applied  for  the  fin  mesh  to  avoid  the  inherent  singularities  of  the  O-topology.  One-half  of  the 
domain  is  meshed  using  Griclgen ,  and  the  other  half  was  generated  using  a  FORTRAN  routine  by 
reflection  along  the  x-z  plane.  In  the  same  manner  discussed  for  the  wing-body  mesh,  the  EPDE 
solver  has  been  run,  applying  the  same  control  functions  on  the  interior  and  the  fin/body  surfaces. 
A  base  line  grid  of  dimensions  43x36x34,  in  the  axial,  normal,  and  bi-normal  directions 
respectively  is  generated. 


*Cr  is  the  Root  Whig- Chord  Length. 

Figure  4-4.  Outer  domain  dimensions 

4.2.3  The  Composite  Mesh.  In  this  subsection  PEGSUS  [51],  the  code  which  implements 
the  Chimera  (Overset)  technique  [49],  is  briefly  discussed,  then  the  process  of  the  creation  of  the 
composite  mesh  will  be  presented.  Figure  4-6  depicts  a  block  diagram  for  the  grid  generation 
process,  starting  from  the  model  geometry  and  ending  up  with  the  composite  mesh.  For  more 
information  on  PEGSUS  and  the  nomenclature  used  throughout  the  rest  of  this  subsection,  the 
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reader  may  refer  to  [49,51].  The  utilized  flow  solver  in  this  study  (the  fdl3di-code)  implements 
the  multi-overset  grid  technique,  and  has  been  successfully  utilized  in  other  applications  [38,47] 
and  validated  [36].  Interpolation  of  conservative  variables  was  adopted  as  recommended  in  the 
literature  [53]. 

PEGSUS  is  designed  to  perform  the  following  tasks  for  each  individual  mesh: 

Establish  the  communication  boundaries  (outer,  hole  or  interior  boundaries). 

Search  for  the  interpolating  stencil  points  in  the  donor  or  recipient  meshes  according  to  the 
established  hierarchical  level. 


Figure  4-5.  Section  in  the  vertical  fin  mesh  demonstrates  the  topology. 

Mark  other  special  points  to  be  excluded  from  the  computational  domain  (interior,  or  orphan 
points). 

At  the  end  of  this  process,  each  point  in  the  mesh  will  be  assigned  the  appropriate  status: 

■  Field  point. 

■  Hole  point. 

■  Boundary  (interpolation)  point. 
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■  Orphan  point. 

■  Interpolating  point. 

■  Interior  point. 

The  embedding  process  of  the  fin  mesh  into  the  wing-body  mesh  requires  the  creation  of  a 
hole  in  the  wing-body  mesh  to  exclude  the  fin  solid  volume  from  the  computational  domain  of  the 
wing-body  mesh.  A  hole-creation  boundary  defined  in  the  fin  mesh  cannot  be  completely  closed; 
a  gap  must  exist  at  the  intersection  of  the  fin  surface  and  the  hole-creation  boundary  (figure  4-7). 
As  a  result,  some  points  close  to  the  fuselage  surface  will  be  considered  unblanked  by  the 
PEGSUS  hole  location  algorithm  (see  figure  4-7).  This  problem  can  be  solved  by  utilizing  the  so 
called  Phantom  mesh,  which  is  a  dummy  mesh  -  in  the  sense  that  it  does  not  communicate  with 
the  other  meshes  in  any  way  —  used  only  for  the  hole  creation.  The  hole-creation  boundary, 
shown  in  figure  4-6,  has  been  used  to  define  the  outer  boundary  of  the  phantom  mesh.  However, 
the  lower  face  of  the  phantom  mesh  is  extended  into  the  fuselage  interior.  The  geometry  of  the 
outer  boundaries  of  the  fin  and  phantom  meshes  is  designed  to  help  capture  the  interpolating 
points  in  the  wing-body  O-grid.  Since  the  fin  mesh  is  an  H-H  grid,  the  fin  solid  surface  should  be 
marked  as  an  interior  region  in  the  fin  mesh  to  be  excluded  from  the  computational  domain  with 
boundary  conditions  being  applied  here.  The  hole  in  the  wing-body  mesh  is  designed  in  a  way  to 
be  interpolated  on  five  faces,  leaving  the  lower  face  for  the  application  of  boundary  conditions. 
The  same  procedure  is  done  for  the  fin  mesh  outer  boundary,  where  the  boundary  conditions  are 
applied  on  the  lower  face.  A  double-fringe  is  globally  applied  for  the  two  meshes  for  better 
results  [51],  and  to  provide  the  required  interpolating  stencil  points  for  the  4*-order  derivatives  of 
the  damping  terms  in  the  flow  solver.  An  overlap  of  approximately  five  points  is  kept  in  the  three 
coordinate  directions  to  satisfy  the  double-fringe  requirement.  It  is  recommended  to  pick  the  in¬ 
coordinate  (or  K)  in  the  axial  direction,  since  PEGSUS  is  designed  to  print  the  MAPS  file  at 
constant  ^-planes.  Also,  it  is  strongly  recommended  to  do  domain  inclusion,  either  globally  or 
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locally,  in  order  to  avoid  any  fault  marking  of  orphan  points  in  the  neighborhood  of  mirror  planes 
(the  z=0  plane  in  our  case).  The  user  input  file,  the  input  mesh  format,  and  the  parameter  values 
required  to  run  PEGSUS  4.0  can  be  found  in  Appendix-A.  A  summary  of  the  PEGSUS  output  is 
found  in  Appendix-B,  showing  no  orphan  points  in  the  two  meshes.  Also,  the  summary  of  the 
output  contains  the  parameter  values  required  to  run  the  fdl3di- flow  solver  (Namely  IGRID, 
ILDIM,  and  IBDIM).  A  sample  of  the  output  maps,  showing  the  layout  of  some  cross-sectional 
planes  in  the  two  meshes,  is  found  in  Appendix-B  under  Diagnostic  Maps. 


Figure  4-6.  Block  diagram  illustrates  the  different  phases  of  the  grid  generation  process. 


Figure  4-7.  The  hole-creation  boundary. 
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Chapter-5  Analysis  and  Results 


5.1  Introduction 

The  main  emphasis  of  this  research  is  to  study  the  effect  of  adding  the  sideslip  degree-of- 
freedom  (DoF)  to  the  analysis  of  wing  rock  onset  and  characteristics.  Since  the  actual  wing  rock 
motion  is  a  combination  of  roll,  sideslip,  and  vertical  descent  (much  like  Dutch-Roll  motion,  see 
figure  5-1),  all  three  of  these  DoF  will  be  included  and  these  three  DoF  solutions  will  be 
compared  and  contrasted  with  the  single  DoF  roll-only  solutions.  Since  the  vertical  fin  can 
produce  a  considerable  amount  of  the  sideslip  force,  the  simulations  are  conducted  for  the  wing- 
body  configuration  with  and  without  the  fin  to  study  its  effect.  All  the  simulations  are  conducted 
for  the  same  initial  flight  condition  (Mach  number  of  0.4  and  pressure  altitude  of  2000  meters).  In 
all  simulations,  the  wing  rock  oscillations  are  triggered  from  an  initial  steady  level  flight  due  to 
an  initial  perturbation  in  roll  rate  (p° )  of  0.01  (non-dimensional  value). 


Figure  5-1.  Schematic  drawing  for  the  wing  rock  motion. 
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In  the  remainder  of  this  Chapter,  some  interesting  phenomenon  will  be  discussed  at 
particular  cross-planes.  For  illustration  the  axial  location  of  these  cross-planes  is  depicted  in 
figure  5-2. 


Figure  5-2.  Axial  location  of  particular  cross-planes. 


5.2  Grid  Sensitivity  Analysis 

Three  meshes  of  different  dimensions  (Table  4-1)  were  developed  to  analyze  the  grid 
dimension  effects  on  the  solution.  The  developed  solutions  for  the  three  meshes  are  compared  and 
analyzed  for  a  sample  run.  The  selected  sample  was  chosen  for  the  three  DoF  motion,  where  the 
dynamics  of  the  motion  depends  on  the  computed  aerodynamic  loads  in  the  three  DoF,  to 
maximize  the  differences  in  solutions.  The  simulations  were  done  for  an  angle  of  attack  (a  =  25°), 
slightly  lower  than  the  expected  wing  rock  onset  (a  =  35°).  The  expected  response  at  a  =  25° 
should  dampen  more  rapidly,  and  hence  reduces  the  cost  of  the  simulation  for  the  three  meshes. 


5-2 


Figure  5-3.  Roll  angle  response  for  the  three  meshes  at  a  =  25°  due  to  an  initial  perturbation  in 

roll  rate  of  0.01. 

Figure  5-3  depicts  the  time-history  for  the  roll  angle  in  response  to  an  initial  perturbation  in 
roll  rate  for  the  three  meshes.  The  results  show  a  convergence  toward  the  fine  mesh  solution,  with 
the  solution  of  the  medium  mesh  very  close  (even  identical  in  some  regions)  to  the  fine  mesh 
solution.  This  supports  that  the  fine  mesh  is  approximately  representing  a  grid-independent 
solution.  The  comparison  of  the  solutions  developed  by  the  coarse  and  fine  meshes  shows  similar 
qualitative  behavior  for  the  coarse  mesh  with  a  7%  increase  in  amplitude  on  average,  and  a 
reduction  in  frequency  of  8%  on  average  over  the  fine  mesh.  The  reason  for  these  differences  is 
the  slightly  stronger  leading  edge  vortex,  which  was  observed  for  the  fine  mesh  (figure  5-4)  at  the 
initial  flight  condition  (25°  angle  of  attack  and  zero  roll  angle).  This  stronger  vortex  increases  the 
damping  slightly  in  the  solution  of  the  fine  mesh  over  that  of  the  coarse  one.  Since  this  research  is 
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primarily  concerned  with  the  incremental  effects  due  to  adding  sideslip  and  vertical  DoF,  it  was 
decided  that  the  coarse  mesh  accuracy  was  adequate  for  this  purpose.  Consequently,  to  save 
computational  cost  the  majority  of  the  simulations  were  conducted  on  the  coarsest  mesh. 

5.3  Simulations  of  Wing  Rock  for  the  Wing-Body  Configuration 

The  geometry  of  the  wing-body  configuration  is  demonstrated  in  detail  in  section  4.1,  which 
represents  the  generic  fighter  model  without  a  fin. 

5.3.1  Single  Degree-of -Freedom  in  Roll.  To  begin  the  research,  an  initial  guess  for  the 
onset  angle  of  attack  was  determined  from  a  previously  published  wing  rock  boundary  versus 
aspect  ratio  [27:733].  The  boundary  in  figure  5-5  was  established  from  the  observations  for  a 
delta  wing  undergoing  wing  rock  in  roll  only.  Although  the  model  in  this  case  is  a  blended  wing- 
body  configuration,  the  predicted  onset  based  on  this  curve  should  be  in  the  neighborhood  of  the 
actual  onset  value.  The  model  aspect  ratio  is  1.29,  which  from  figure  5-5  corresponds  to  an  onset 
of  33  degrees  in  angle  of  attack.  Also  the  observations  in  [27]  indicate  that  the  wing  rock  is  likely 
to  be  triggered  and  driven  by  asymmetric  vortex  bursting,  and  not  by  asymmetric  vortex 
shedding.  The  33-degrees  angle  of  attack  was  run  for  a  couple  of  cycles,  and  didn’t  lead  to  a  limit 
cycle  oscillation.  Next  a  35-  degrees  angle  of  attack  was  tried  and  successfully  led  to  a  limit  cycle 
oscillation  in  roll  (figure  5-6).  Approximately  16  degrees  peak-to-peak  amplitude  in  roll  angle  is 
observed  with  an  interesting  limit  cycle  phase-plane  portrait  (figure  5-7). 
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Figure  5-7.  The  phase-plane  portrait  for  a  full  limit  cycle  in  roll. 

The  relatively  low-amplitude  of  the  wing  rock  oscillation  was  actually  not  surprising,  the 
previous  work  of  Brandon  and  Nguyen  [4:3]  on  a  chine  fore-body  mounted  on  a  26°-sweep 
trapezoidal  wing  found  oscillations  of  approximately  12°  amplitude  (see  figure  1-3).  A  rapid 
rolling  rate  was  observed  in  the  neighborhood  of  the  zero-roll  angle  position  (look  at  the  little 
circled  kink  in  figure  5-6).  Surprisingly,  this  little  kink  in  the  roll  angle  response  has  been 
observed  in  the  flight  test  data  of  wing  rock  for  realistic  fighters.  The  closest  configuration  to  the 
generic  fighter  model  geometry,  for  which  wing  rock  flight  test  data  is  available,  is  the  F-4  fighter 
aircraft  [23]  (figure  5-8).  The  recorded  flight  test  data  for  the  F-4  wing  rock  [43:6-13],  shows  a 
similar  roll  angle  response  (figure  5-9)  to  that  of  the  generic  fighter  model  (figure  5-6). 
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Roll  Angie  (degrees) 


Figure  5-9  Roll  angle  response  following  the  onset  of  wing  rock  for  the  F-4  [43:6-13] 


In  order  to  explain  the  fluid  mechanism  that  drives  the  observed  wing-body  rock  in  this 
study,  two  main  phenomena,  the  vortex  burst  and  vortex  dynamics,  should  be  investigated 
throughout  a  full  cycle  of  wing  rock. 

The  instantaneous  axial  location  of  vortex  burst  was  traced  by  searching  the  first  occurrence 
of  a  negative  axial  velocity  on  the  upper  surface  of  the  two  wing  halves.  The  instantaneous  axial 
locations  of  vortex  burst  captured  by  the  searching  algorithm  within  an  entire  cycle  of  wing  rock 
are  depicted  in  figure  5-10.  In  order  to  verify  the  instantaneous  axial  locations  of  vortex  burst 
captured  by  the  searching  algorithm,  mass-less  particles  were  injected  near  the  tips  of  the  leeward 
side  of  the  body.  The  traced  trajectories  of  the  particles  at  the  initial  conditions  ( a  =  35°,  cp  =  0°) 
showed  a  flow  reversal  (see  figure  5-11)  at  approximately  the  same  predicted  location  from  the 
numerical  algorithm  (x*=  1.032). 

In  order  to  explain  the  vortex  burst  dynamics  depicted  in  figure  5-10,  let  us  assume  the 
model  is  rolling  in  the  positive  direction.  Then  at  a  positive  roll  angle,  the  axial  location  of  vortex 
burst  on  the  left  wing  (which  moves  upward)  moves  downstream,  while  that  on  the  right  wing 
(which  moves  downward)  is  moving  upstream.  In  the  light  of  this  behavior  of  the  asymmetric 
vortex  burst  location,  one  can  explain  its  role  in  driving  the  oscillations.  Once  the  model  is 
initially  perturbed  in  roll  rate,  it  starts  to  oscillate  around  the  zero-roll  angle  position.  The 
upward-moving  wing  is  exposed  to  a  down  wash,  which  reduces  the  instantaneous  value  of  angle 
of  attack  on  this  wing  and  hence  delays  vortex  burst  (i.e.  the  first  occurrence  of  vortex  burst 
moves  downstream).  On  the  other  hand,  the  downward-moving  wing  is  exposed  to  an  up  wash, 
which  increases  the  instantaneous  value  of  angle  of  attack  on  this  wing,  and  hence  promotes 
vortex  burst  on  this  wing  (i.e.  the  first  occurrence  of  vortex  burst  moves  upstream).  As  a  result, 
the  lift  is  increased  on  the  upward-moving  wing  and  reduced  on  the  downward-moving  wing.  So 
far,  this  model  should  go  into  a  divergent  rolling  motion  in  the  absence  of  a  restoring  moment. 
However,  the  effect  of  the  variation  of  vortex  structure  and  location  (vortex  dynamics)  still  need 
to  be  investigated  to  reach  a  complete  understanding  of  the  physics  of  the  motion.  The  vortex 
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dynamics  was  visualized  in  this  study  through  snap  shots  of  the  instantaneous  stagnation  pressure 
coefficient  distribution  for  a  cross-plane  at  x*  =  1.06  (see  figure  5-12).  As  shown  in  figure  5-12, 
the  vortex  on  the  upward-moving  wing  becomes  more  coherent  and  moves  slightly  away  of  the 
surface,  so  its  footprint  on  the  surface  is  reduced  (see  figure  5-13  (b)).  While  the  vortex  on  the 
downward  moving  wing  is  diffused  and  slightly  moves  toward  the  surface,  so  its  footprint  on  the 
surface  is  increased  (see  figure  5-13  (b)).  As  a  result,  the  lift  is  reduced  on  the  upward-moving 
wing  and  increased  on  the  downward-moving  wing,  which  establishes  the  restoring  moment. 
Figure  5-13  (a)  shows  the  symmetric  suction  on  the  wing  upper  surface  at  the  initial  conditions  (a 
=  35°,  (p  =  0°),  while  figure  5-13  (b)  shows  an  asymmetric  suction  on  the  wing  upper  surface  at  a 
roll  angle  of  8°  due  to  vortex  dynamics.  The  effect  of  vortex  dynamics  during  wing-body  rock  is 
counter-balancing  the  effect  of  the  asymmetric  vortex  burst  location,  which  sustains  the  wing 
rock  oscillation. 

In  summary,  the  vortex  burst  dynamics  is  found  to  produce  the  driving  moment  of  the 
motion  for  small  roll  angles,  where  the  asymmetry  in  vortex  structure  and  location  is  not 
significant.  As  the  rolling  continues,  the  rolling  angle  increases  and  the  asymmetries  in  vortex 
structure  and  location  becomes  significant,  and  hence  the  restoring  moment  increases.  At  a  finite 
roll  angle,  the  restoring  moment  becomes  sufficiently  strong  to  stop  rolling  and  even  reverses  the 
direction  of  roll  afterwards. 

Based  on  the  previous  explanation  of  the  fluid  mechanism  that  drives  wing  rock,  tangential 
blowing  on  the  upper  surface  to  enhance  vortex  burst  is  suggested  for  wing  rock  delay  or 
suppression.  Forebody  vortex  control  has  been  studied  for  a  similar  model  consisting  of:  forebody 
of  circular  cross-section,  sharp  delta  wing  of  78°-sweep,  and  vertical  fin  [41].  The  reported  results 
showed  that  steady  tangential  blowing  from  leeward  nozzles  near  the  forebody  tips  was  capable 
of  suppressing  wing  rock  [41:303]. 
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The  previous  discussion  explains  the  interesting  shape  of  the  limit  cycle  (figure  5-7)  due  to 
the  rapid  increase  in  roll  rate  in  the  neighborhood  of  the  zero-roll  angle  position.  Since  the 
asymmetries  in  vortex  structure  and  location  in  the  neighborhood  of  the  zero-roll  angle  position 
are  not  significant,  the  asymmetric  vortex  burst  location  is  dominantly  driving  the  motion.  This 
also  explains  the  switch  in  the  sign  of  the  roll  stiffness  derivative  C/^from  negative  to  positive  as 
the  model  passes  points  (A)  and  (B)  (see  figure  5-14). 


Roll  Angle  [deg] 


O  Negative  Roll  Rate 

+  Positive  Roll  Rate 

- - Wings-Leveled  Static  Location 

- Left  Wing” 

-  Right  Wing" 


"The  Left  &  Right  Conventions  are  Referred  to  an 
Observer  Standing  at  the  Origin  in  the  Y-Z  Plane 
and  Looking  down  the  (+)  direction  of  X-axis. 


Figure  5-10.  Instantaneous  vortex  burst  locations  during  an  entire  cycle  of 
roll  oscillations  at  a  =  35°. 

Figure  5-10  shows  large  hysteresis  associated  with  the  vortex  burst  location  on  both  wing 
halves.  This  hysteresis  is  observed  on  the  wing-half  only  if  it  is  in  the  leeward  side  of  the  wind. 
Even  though  the  model  discussed  here  does  not  have  sideslipping  body  motion,  there  is  still  a 
kinematically  induced  sideslip.  Since  the  model  is  rolling  about  the  body-fixed  x-axis,  the  angle 
of  attack  and  sideslip  are  mutually  interchanging  during  rolling  as  follows: 
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tan  p  =  tan  an  sin  (p 


(5-1) 


Where: 

P  =  Sideslip  angle  due  to  kinematic-coupling. 
a„  =  Angle  of  attack  (wings-leveled). 

<p  =  Roll  angle. 

At  a  positive  roll  angle  the  sideslip  is  positive  and  hence  the  left  wing  is  in  the  leeward  side 
of  wind.  Similarly  at  a  negative  roll  angle,  the  sideslip  is  negative  and  hence  the  right  wing  is  in 
the  leeward  side  of  wind.  The  hysteresis  associated  with  the  instantaneous  vortex  burst  location 
(figure  5-10)  is  only  observed  on  the  upward-moving  wing  due  to  the  time-lag  of  the  vortex 
motion. 

5.3.2  Three  Degrees-of  -Freedom  Wing  Rock.  The  previous  simulation  of  wing  rock  for 
the  wing-body  configuration  with  a  single  DoF  in  roll  at  a  =  35°  is  next  repeated  for  the  three 
DoF:  roll,  sideslip,  and  vertical  descent.  The  simulation  shows  the  response  starts  divergent,  and 
damps  out  nicely  after  the  first  cycle  (see  figure  5-15).  This  result  indicates  that  adding  the 
sideslip  DoF  increases  the  damping  of  the  motion  in  general  and  hence  delays  the  onset. 
However,  the  simulation  should  be  conducted  to  a  higher  value  of  angle  of  attack  to  confirm 
whether  the  system  would  go  into  a  limit  cycle  or  not  at  all.  Raising  the  angle  of  attack  above  the 
onset  (a  =  35°)  one  degree  at  a  time  is  computationally  costly  (at  least  3-4  cycles  are  needed  to 
capture  a  limit  cycle),  so  the  angle  of  attack  was  raised  by  5°.  The  simulation  was  repeated  for  a 
=  40°  for  the  same  initial  conditions  (p°  =  .01),  and  now  achieves  a  limit  cycle  (figure  5-16).  This 
result  concludes  that  adding  the  sideslip  DoF  delays  the  onset  of  wing  rock  in  general,  however 
this  doesn’t  mean  that  the  onset  was  delayed  by  5°,  since  the  angles  between  35°  and  40°  have  not 
been  tested.  Comparing  figures  5-6  and  5-16  one  can  see  that  the  amplitude  of  the  oscillations 
was  reduced  by  approximately  50%  over  the  single  DoF  case. 
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Figure  5-11  Particle  path  trajectories  nt  35  degree 


(d)  cp  =  -3.5°  (h)  (p  =  3.0° 

Figure  5-12.  Vortex  dynamics  during  a  full  limit  cycle  in  roll  at  a  =  35°. 
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Figure  5-1 


(b)  <p  =  8° 


3.  Pressure  coefficient  distribution  at  a  =  35°  for  a  cross-section  at  /=!  .06. 
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Figure  5-14.  Variations  of  rolling  moment  coefficient  versus  roll  angle 
for  a  full  limit  cycle  in  roll. 

An  aircraft  in  a  rolling  motion  about  its  body-fixed  longitudinal  axis  is  subjected  to  sideslip 
due  to  two  different  sources,  one  from  the  kinematic-coupling  of  the  angle  of  attack  and  sideslip 
(equation  5-1)  and  the  other  from  the  dynamics  of  the  body’s  sideward  motion  (i.e.  kinetics).  The 
sideslip  due  to  the  kinematic-coupling  is  in  phase  with  the  motion  as  shown  in  figure  5-17,  while 
the  sideslip  due  to  kinetics  lags  the  motion.  Figure  5-17  shows  an  approximately  90°  phase-lag 
between  the  rolling  and  sideslip  due  to  kinetics.  This  90°  phase-lag  is  a  source  of  adverse  sideslip 
motion  during  rolling,  which  has  the  effect  of  reducing  the  vortex  lift-off  distance  on  the  upward- 
moving  wing  and  hence  reducing  the  amplitude  of  the  oscillations.  In  addition,  this  adverse 
sideslip  has  the  effect  of  reducing  the  hysteresis  associated  with  the  instantaneous  vortex  burst 
location  (figure  5-18)  compared  with  the  case  for  the  single  DoF  (figure  5-10). 
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Nondimensional  Time 


Nondimensional  Time 


Figure  5-17.  Roll  and  sideslip  angles  time-histories  for  an  entire  cycle  of  wing  rock 
at  a  =  40°  for  the  case  of  three  DoF  motion. 


Roll  Angle  [deg] 


Figure  5-18.  Instantaneous  vortex  burst  locations  during  an  entire 
cycle  of  the  three  DoF  oscillations  at  a  =  40°. 


Previously,  it  has  been  shown  that  the  static  lateral  stability  derivative  (known  as  the 

dihedral  derivative)  can  significantly  influence  the  onset  of  wing  rock  [33,34,35,48].  The  dihedral 
derivative  (as  are  all  stability  derivatives)  is  determined  at  a  certain  angle  of  attack  by  measuring 
the  value  of  the  rolling  moment  coefficient  for  several  static  sideslip  angles  and  then  computing 
the  derivative  from  this  data.  Due  to  computational  costs,  only  dynamic  simulations  were 
performed  for  the  present  research.  However,  one  can  look  at  a  pseudo  stability  derivative  for  the 
dynamic  case  from: 

r 

dp 

dC , 

_  dt 
d_l 
dt 

Where: 

Ci  =  Rolling  moment  coefficient. 

P  =  Sideslip  angle  (sum  of  kinematic-coupling  plus  kinetic  motion)  [rad]. 

The  derivative  was  numerically  computed  by  a  second  order  accurate  central  difference 
scheme  from  the  recorded  data  for  a  full  cycle  of  wing  rock  in  three  DoF.  The  computed  C,  was 


then  filtered  through  a  low-pass  filter  to  suppress  the  high  frequency  noise  associated  with  the 
numerical  scheme.  Figure  5-19  shows  the  raw  and  filtered  data.  The  mean  value  of  Clp  over  a  full 

cycle  was  calculated  to  be  -0.0591  [rad1].  This  concludes  that  during  a  fully  developed  wing  rock 
oscillations,  the  mean  value  of  the  dihedral  derivative  is  strongly  stable.  This  contradicts  the 
results  of  reference  [48],  where  it  was  shown  for  the  F-15  that  C,  is  approximately  zero  at  the 

onset  of  wing  rock.  It  remains  to  be  seen  whether  this  pseudo  stability  derivative  is  an  accurate 
reflection  of  the  true  static  definition.  Further  numerical  static  simulations  would  have  to  be 
performed. 
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5.4  Simulations  of  Wing  Rock  for  the  Full  Generic  Fighter  Model 

The  fin  can  produce  a  considerable  amount  of  the  side  force,  so  its  effect  on  the  three  DoF 
simulations  needs  to  be  examined.  A  finite-thickness  trapezoidal  fin  of  moderate  fin-to-wing  area 
ratio  of  0.08  was  used.  The  simulation  of  wing  rock  at  a  =  35°  was  repeated  for  the  generic 
fighter  model  including  the  fin.  The  results  show  an  increase  in  the  damping  and  frequency  of  the 
oscillations  (figure  5-20).  However,  this  result  is  limited  to  the  fin  geometry  used  in  this  study, 
since  the  effect  of  a  fin  is  greatly  dependent  upon  its  geometry  and  its  area  ratio  [37:500].  The 
analysis  reported  in  [41:301,4:3],  for  the  case  of  single  DoF  in  roll,  shows  that  fins  of  relatively 
small  area  ratios  have  the  effect  of  delaying  the  onset  and  reducing  the  amplitude  of  wing  rock. 
Also  it  has  been  reported  that  fins  of  large  area  ratios  are  likely  to  get  separated  and  produce  an 
opposite  effect  [41:301].  The  three  DoF  simulation  of  wing  rock  including  the  fin  shows  a 
significant  increase  in  frequency  (figure  5-20).  This  has  not  been  reported  for  most  of  the  fin 
configurations  in  literature  for  the  case  of  single  DoF  roll  only  wing  rock.  Only  one  double  fin 
configuration  has  shown  a  similar  effect  for  the  single  DoF  in  roll  [37:500].  A  trapezoidal  plane- 
form  of  area  ratio  .005  led  to  an  increase  in  frequency  of  36%,  and  significant  reduction  in  wing 
rock  amplitude  [37:500]. 

The  fin  has  two  main  contributions,  one  due  to  the  aerodynamic  interaction  between  the  fin 
surface  and  the  flow  over  the  wing-body,  and  the  other  due  to  the  side  force  produced  by  the  fin. 
The  main  observed  aerodynamic  effect  of  the  fin  was  the  significant  delay  of  vortex  burst.  The 
instantaneous  location  of  vortex  burst  was  determined  using  the  same  criterion  that  was 
previously  explained  in  subsection  5.3.1.  Figure  5-21  shows  that  the  vortex  burst  location  at  a  = 
35°  for  the  model  with  fin  was  first  observed  at  x*  =  1.48,  just  downstream  of  the  fin  trailing  edge 
compared  with  x*  =  1.03  for  the  model  without  fin.  Figure  5-22  shows  the  non-dimensional  axial 
velocity  distribution  at  x  =  1.44  (slightly  before  the  fin  trailing  edge),  which  shows  a  positive 
minimum  axial  velocity.  So  vortex  bursting  has  not  been  observed  upstream  of  the  trailing  edge 
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of  the  fin  at  a  =  35°.  Since  vortex  burst  dynamics  produces  the  driving  moment  for  the 
oscillations,  the  oscillations  damped  out  more  rapidly  with  the  addition  of  the  fin. 


Figure  5-20.  Roll  angle  response  to  an  initial  roll  rate  of  0.01  at  a  =  35' 
in  three  DoF  for  the  model  with  and  without  vertical  fin. 


Figure  5.21  Cross-plane  axial  velocity  distribution  at*  =  1.48 
for  the  generic  fighter  at  a  =  35°. 


Figure  5.22  Cross-plane  axial  velocity  distribution  at  x'  =  1 .44 
for  the  generic  fighter  at  a  =  35°. 


5.5  Computational  Cost 

The  CPU- time  for  the  diagonalized-Euler  version  of  the  modified  fdl3di-c ode  for  the  case  of 
three  DoF  is  approximately  2.51  /rs/grid-point/sub-iteration  on  the  CRAY-T94  supercomputer. 
Negligible  difference  was  observed  between  the  computational  cost  for  the  case  of  single  DoF  in 
roll  and  that  for  the  case  of  three  DoF.  The  reason  is  that  the  computational  cost  of  solving  the 
rigid-body  dynamics  equations  is  negligible  compared  with  the  cost  of  solving  the  flow  governing 
equations.  This  fact  motivates  conducting  the  computational  simulations  of  wing  rock  in  three 
DoF  rather  than  single  DoF  for  better  accuracy. 
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Chapter-6  Conclusions  and  Recommendations 


In  the  present  study,  the  lateral  limit  cycle  oscillation  phenomenon  known  as  wing  rock  has 
been  analyzed  in  single  and  three  degrees-of-freedom  (DoF)  for  a  generic  fighter  model  with  and 
without  a  vertical  fin. 

The  results  for  the  model  without  fin  in  single  DoF  in  roll  showed  a  wing  rock  onset  at  a  = 
35°  with  several  interesting  phenomena.  A  roll  amplitude  of  16°  (peak-to-peak)  was  observed, 
which  confirms  the  early  results  [4]  on  a  similar  model.  The  dynamic-aerodynamic  interactions 
during  the  motion  were  thoroughly  investigated  by  analyzing  the  snap  shots  taken  for  the 
flowfield  at  specific  time  intervals  and  tracing  the  instantaneous  location  of  vortex  burst  on  both 
wing  halves  for  an  entire  cycle.  The  vortex  burst  dynamics  during  the  motion  was  found  to 
produce  the  driving  moment,  while  the  vortex  dynamics  (lift-off)  near  the  peaks  of  the 
oscillations  on  the  upward-moving  wing  half  restores  the  motion. 

For  comparison  with  the  single  DoF  in  roll  simulations,  simulations  were  conducted  in  three 
DoF:  roll,  sideslip,  and  vertical  motion.  The  purpose  of  the  three  DoF  simulations  was  to  study 
the  effect  of  the  sideslip  dynamics  on  the  onset  and  limit  cycle  characteristics  of  wing  rock.  The 
results  on  the  model  without  fin  demonstrated  a  delay  in  the  onset  and  a  high  reduction 
(approximately  50%)  in  the  amplitude  of  the  oscillations  with  surprisingly  no  change  in 
frequency.  The  associated  sideslip  kinetics  was  found  to  lag  the  rolling  motion  by  90°  during  the 
fully  developed  wing  rock  oscillations.  This  adverse  sideslip  motion  was  found  to  be  the  source 
of  this  augmentation  in  damping  and  delay  in  onset.  Also,  the  mean  value  of  a  pseudo  dihedral 
derivative  during  a  full  cycle  of  oscillation  was  found  to  be  clearly  stable  (-.0591),  which 
contradicts  previously  reported  F-15  results  [48]. 

The  effect  of  the  vertical  fin  on  the  onset  and  characteristics  of  wing  rock  in  three  degrees- 
of-freedom  was  found  similar  to  the  results  reported  in  [4,37,41]  for  the  case  of  single  degree-of- 
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freedom  in  roll.  The  vertical  fin  was  observed  to  greatly  affect  the  frequency  for  the  case  of  three 
DoF  wing  rock,  while  no  appreciable  change  in  frequency  was  observed  for  the  case  of  single 
DoF  roll  only  wing  rock.  The  present  research  shows  that  the  vertical  fin  enhances  the  vortex 
burst  on  the  wing  by  delaying  its  location  downstream  of  the  fin  trailing  edge  and  hence  preserves 
the  flow  energy  on  the  wing. 

In  summary,  the  following  conclusions  can  be  drawn  concerning  wing  rock  for  the  generic 
fighter  model  described  earlier: 

The  analysis  of  the  wing  rock  motion  for  a  full  cycle  showed  that  the  vortex  burst 
dynamics  is  the  source  of  the  driving  moment  that  sustains  the  motion,  while  the  vortex 
lift-off  on  the  up-ward  moving  wing  creates  the  restoring  moment. 

In  computational  simulations,  the  sideward  and  vertical  DoF  should  be  included  to  get  an 
accurate  prediction  of  the  wing  rock  onset  and  limit  cycle  characteristics.  Especially,  as 
previously  mentioned  in  section  5.5,  the  computational  cost  of  including  the  other  two 
DoF  is  negligible.  However,  the  predicted  values  for  the  single  DoF  in  roll  only  were 
found  to  be  conservative  (i.e.  earlier  onset  prediction  and  larger  amplitude  limit  cycle 
characteristics). 

The  computed  value  of  the  mean  pseudo  dihedral  derivative  for  a  full  cycle  of  the  fully 
developed  three  DoF  simulations  of  wing  rock  was  found  clearly  stable. 

-  Discarding  the  vertical  fin  in  the  prediction  of  the  wing  rock  onset  and  characteristics  was 
found  to  lead  to  different  results.  Unlike  the  effect  of  including  the  three  DoF  in  wing  rock 
simulations,  which  is  approximately  model-independent,  the  effect  of  vertical  fin  is 
greatly  depending  on  its  geometry  and  can’t  be  generalized. 
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The  author  recommends  performing  the  following  utilizing  the  CFD-code  developed  for  this 
effort: 

-  Computing  the  static  lateral  stability  derivatives  at  the  predicted  onset  to  correlate  with  the 
pseudo  derivative  that  was  computed  from  the  recorded  data  during  the  motion. 

Studying  the  effect  of  fore-body  vortex  control  by  tangential  blowing  for  the  case  of  three 
DoF  wing  rock  oscillations. 

-  Obtain  more  accurate  results  for  the  onset  point  by  repeating  the  simulations  for  1° 
increments  in  angle  of  attack  rather  than  the  5°  increment  used  in  the  present  study. 

Perform  more  fine  mesh  runs  to  validate  the  coarse  mesh  results. 

Extending  the  present  work  for  flexible-body  dynamics. 
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Appendix  A  -  PEGSUS  Inputs 


A.l  PEGSUS  User  Inputs  File 

PEGSUS  4.0  INPUT  FILE  (Coarse  Mesh  Version) 

DEFINE  GLOBAL  PARAMETERS 
$ GLOBAL 
FRINGE  =  2, 

$END 

C  DEFINE  MESHES 

$MESH  NAME  =  ' CWBGRID 1 , 

LINK  =  ' CFINGRID 1 , 

JINCLUDE  =  31,68, 

KINCLUDE  =  2,28, 

LINCLUDE  =  19,47, 

QUALITY  =  1.0, 0.3, -0.1, 

$END 

$MESH  NAME  =  ' CFINGRID ' , 

LINK  =  ' CWBGRID ' , 

JINCLUDE  =  1,43, 

KINCLUDE  =  1,36, 

LINCLUDE  =  1,34, 

QUALITY  =  1.0, 0.3, -0.1, 

$END 

$MESH  NAME  =  ' PHANTGRID ' 

$END 

C  DEFINE  BOUNDARIES 

$ BOUNDARY  NAME  =  ' PHANT  HOLE  BOUND ' , 

ISPARTOF  =  'PHANTGRID', 

MHOLEIN  =  ' CWBGRID ' , 

$END 

$BOUNDARY  NAME  =  ' FIN  OUTER  BOUND ' , 

ISPARTOF  =  'CFINGRID', 

CLOSED  =. FALSE., 

$END 

C  DEFINE  THE  SURFACES  THAT  CONNECTED  TO  THE  PHANT  HOLE  BOUND 

C  THE  BOUNDARY  IN  GENERAL  IS  A  BOX- LIKE  SHAPE  COMPRISES  OF  4 

C  SURFACES  TO  AVOID  CONCAVITIES  AT  THE  BOX  CORNERS 

$  SURFACE  ISPARTOF  =  'PHANT  HOLE  BOUND', 

J RANGE  =  1,40, 

KRANGE  =  1,10, 
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L RANGE  =  30,30, 
NVOUT  =  ' +L '  , 
$END 

$ SURFACE  ISPARTOF 
J RANGE  =  1,40, 
KRANGE  =  1,10, 

L RANGE  =  1,1, 
NVOUT  =  ' -L ' , 
$END 

$ SURFACE  ISPARTOF 
J RANGE  =  1,1, 
KRANGE  =  1,10, 

L RANGE  =  1,30, 
NVOUT  =  ' -J ' , 
$END 

$ SURFACE  ISPARTOF 
J RANGE  =  40,40, 
KRANGE  =  1,10, 

L RANGE  =  1,30, 
NVOUT  =  ' + J ' , 
$END 

$ SURFACE  ISPARTOF 


JRANGE  = 

1,40, 

KRANGE  = 

1,1, 

LRANGE  = 

1,30, 

NVOUT  = 

’-K'  , 

$END 

$ SURFACE 

ISPARTOF 

JRANGE  — 

1,40, 

KRANGE  = 

10,10, 

LRANGE  = 

1,30, 

NVOUT  = 

'  +K '  , 

$END 


'  PHANT  HOLE  BOUND ' , 


'  PHANT  HOLE  BOUND ' , 


1  PHANT  HOLE  BOUND '  , 


1  PHANT  HOLE  BOUND ' , 


'  PHANT  HOLE  BOUND ' , 


DEFINE  THE  OUTER  BOUNDARY  OF  THE  FIN  MESH 


$ SURFACE  ISPARTOF 
JRANGE  =  1,43, 
KRANGE  =  1,36, 

L RANGE  =  33,34, 
$END 

$ SURFACE  ISPARTOF 
JRANGE  =  1,2, 
KRANGE  =  1,36, 
LRANGE  =  2,32, 
$END 

$ SURFACE  ISPARTOF 
JRANGE  =  42,43, 
KRANGE  =  1,36, 
LRANGE  =  2,32, 
$END 


'  FIN  OUTER  BOUND ' , 


'  FIN  OUTER  BOUND  1  , 


'  FIN  OUTER  BOUND  1 , 


'  FIN  OUTER  BOUND '  , 


$ SURFACE  ISPARTOF  = 

J RANGE  =  3,41, 

KRANGE  =  1,2, 

L RANGE  =  2,32, 

$END 

$ SURFACE  ISPARTOF  =  ’FIN  OUTER  BOUND’, 

J RANGE  =  3,41, 

KRANGE  =  35,36, 

L RANGE  =  2,32, 

$END 

C  Define  the  Fin  surface  and  the  fuselage  part  in  the  Fin  mesh 

C  as  interior  banked  points  to  be  excluded  from  the  flow  field 
$ REG I ON  NAME  =  'FIN  REGION', 

TYPE  =  ' INTR ' , 

ISPARTOF  =  ' CFINGRID ' , 

$END 

C  DEFINE  FIN  SURFACE 

$VOLUME  ISPARTOF  =  ' FIN  REGION ’ , 

J RANGE  =  14,29, 

KRANGE  =  18,19, 

L RANGE  =  2,20, 

$END 


A. 2  PEGSUS  Parameters 

MDIM=3 

LNDIM=3 

ICHAR=40 

NBDIM=3 

MHDIM=2 

NSDIM=14 

MSLEN=25000 

NXDIM=1 

NRDIM=2 

NVDIM=3 

INCORE=0 

MLEN=999999 

MLEMAX=  999999 

MILEN=250000 

NBYTEI=8 

NBYTER=  8 

MSSMAX  =  MSLEN* (INCORE*NSDIM+l-INCORE) , 
MIMAX1  =  MILEN* (INCORE*MDIM+l-INCORE) , 
MIMAX2  =  MILEN* ( INCORE *MDIM*MDIM+1- INCORE) 
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A. 3  Input  Mesh  Format 

The  individual  meshes  are  input  to  the  PEGSUS  via  the  input 
file  INGRID,  which  has  the  following  format  for  each  individual 
mesh : 

RECORD  VARIABLES 

1  Mesh  Name 

2  JMAX,  KMAX,  LMAX 

3  ( ( (X ( J, K, L) , J=l, JMAX) , K=1 , KMAX) ,L=1,LM 

AX)  , 

( ( (Y ( J, K, L) , J=l, JMAX) , K=1 , KMAX) ,L=1,LM 
AX)  , 

( ( ( Z ( J, K, L) , J=l, JMAX) , K= 1 , KMAX ) ,L=1,LM 
AX) 

Where  JMAX,  KMAX,  and  LMAX  are  the  maximum  mesh  indices  in  the  J, 
K,  and  L  directions  respectively. 
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Appendix  B  -PEGS US  Output 


B.l  PEGS  US  Output  Summary 

PPPPPPPPP  EEEEEEEEE  GGGGGGGGG  SSSSSSSSS  UUU  UUU  SSSSSSSSS 

PPPPPPPPP  EEEEEEEEE  GGGGGGGGG  SSSSSSSSS  UUU  UUU  SSSSSSSSS 
PPP  PPP  EEE  GGG  SSS  UUU  UUU  SSS 

PPP  PPP  EEE  GGG  SSS  UUU  UUU  SSS 

PPPPPPPPP  EEEEEEEEE  GGG  SSSSSSSSS  UUU  UUU  SSSSSSSSS 

PPPPPPPPP  EEEEEEEEE  GGG  SSSSSSSSS  UUU  UUU  SSSSSSSSS 

ppp  EEE  GGG  GGGG  SSS  UUU  UUU  SSS 

PPP  EEE  GGG  GG  SSS  UUU  UUU  SSS 

ppp  EEEEEEEEE  GGGGGGGGG  SSSSSSSSS  UUUUUUUUU  SSSSSSSSS 

PPP  EEEEEEEEE  GGGGGGGGG  SSSSSSSSS  UUUUUUUUU  SSSSSSSSS 


VERSION  4.x 


PARAMETERS 

MDIM=  3 
LNDIM=  3 

ICHAR=  40 
NBDIM=  3 

MHDIM=  2 

NSDIM=  14 
MSLEN=  25000 
NXDIM=  1 

NRDIM=  2 

NVDIM=  3 

INCORE=  0 
MLEN=  999999 
MLEMAX=  999999 

MILEN=  250000 


GLOBAL  DESCRIPTION  READ 
3  MESH  DESCRIPTIONS  READ 
2  BOUNDARY  DESCRIPTIONS  READ 
11  SURFACE  DESCRIPTIONS  READ 
0  BOX  DESCRIPTIONS  READ 
1  REGION  DESCRIPTIONS  READ 
1  VOLUME  DESCRIPTIONS  READ 


READING  MESH:  CWBGRID 
READING  MESH:  CFINGRID 
READING  MESH:  PHANTGRID 
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*****  GLOBAL  PARAMETER  SUMMARY 


*  *  *  *  * 


MIN/MAX  OF  DOMAIN 

MIN  X  =  -0 . 1200E+01 

MAX  X  =  0 . 3000E+01 

MIN  Y  =  -0 . 1500E+01 

MAX  Y  =  0 . 1500E+01 

MIN  Z  =  -0 . 1500E+01 

MAX  Z  =  0 . 1500E+01 

NO  INCLUDE  RANGES  SPECIFIED 

NO.  OF  FRINGE  POINTS  = 

2 

QUALITY { 1 ) 

1.0000 

QUALITY ( 2 ) 

1.0000 

QUALITY ( 3 ) 

1.0000 

EPS 

0 . 1000E-02 

GLOBAL  ROTATION  SEQUENCE  = 

1  *****  MESH  PARAMETER  SUMMARY  ***** 


MESH  NAME:  CWBGRID 

PEGSUS  MESH  NO.=  1 
COMPOSITE  MESH  NO . =  1 


INDEX  RANGES 

:  JMAX  = 

79  KMAX 

=  35 

LMAX  =  65 

TRANSLATION : 

XO 

0 . 0000E+00 

Y0 

0 . 0000E+00 

ZO 

ROTATION: 

XR 

0 . 0000E+00 

YR 

0 . 0000E+00 

ZR 

SCALING: 

ALPHA  = 

SCALE  = 

0 . 0000E+00 

0 . 1000E+01 

BETA  = 

0 . 0000E+00 

GAM 

ROTATION  SEQUENCE  = 

NO.  OF  FRINGE  POINTS 
QUALITY ( 1 ) 

QUALITY (2) 

QUALITY (3) 

EPS 


1.0000 
0.3000 
-0.1000 
0 . 1000E-02 


MIN/MAX  OF  MESH 


MIN 

X  = 

-0 . 1200E+01 

MAX  X  = 

0 .3000E+01 

MIN 

Y  = 

“0 . 1500E+01 

MAX  Y  = 

0 . 1500E+01 

MIN 

Z  = 

-0.1500E+01 

MAX  Z  = 

0 . 1500E+01 

INCLUDE  PARAMETERS 


MIN 

X  = 

-0 . 1000E+31 

MAX  X  = 

0 .1000E+31 

MIN 

Y  = 

-0 . 1000E+31 

MAX  Y  = 

0 .1000E+31 

MIN 

Z  = 

-0 . 1000E+31 

MAX  Z  = 

0 . 1000E+31 

MIN 

J  = 

31 

MAX  J  = 

68 

MIN 

K  = 

2 

MAX  K  = 

28 

MIN 

L  = 

19 

MAX  L  = 

47 
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0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 


LINKS 

NO .  NAME 

1  CFINGRID 


MESH  NAME:  CFINGRID 


PEGSUS  MESH  NO.=  2 
COMPOSITE  MESH  NO . =  2 


INDEX  RANGES: 

:  JMAX  = 

43  KMAX 

=  36 

LMAX  =  34 

TRANSLATION: 

ROTATION: 

XO 

XR 

ALPHA  = 

0 . 0000E+00 
0. 0000E+00 

0 . 0000E+00 

Y0 

YR 

BETA  = 

0 .OOOOE+OO 

0 . 0000E+00 

0. 0000E+00 

zo  = 

ZR  = 

GAM  = 

0. 0000E+00 
0 .0000E+00 
0. 0000E+00 

SCALING:  SCALE  =  0.1000E+01 

ROTATION  SEQUENCE  = 


NO.  OF  FRINGE  POINTS 
QUALITY ( 1 ) 

QUALITY ( 2 ) 

QUALITY (3) 

EPS 

MIN/MAX  OF  MESH 
MIN  X  =  0 . 9400E+00  MAX  X  =  0.1720E+01 

MIN  Y  =  -0 . 4000E+00  MAX  Y  =  0.4000E+00 

MIN  Z  =  0 . OOOOE+OO  MAX  Z  =  0.4599E+00 

INCLUDE  PARAMETERS 
MIN  X  =  -0 . 1000E+31 

MIN  Y  =  -0 . 1000E+31 

MIN  Z  =  -0 . 1000E+31 

MIN  J  =  1 

MIN  K  =  1 

MIN  L  =  1 


LINKS 

NO .  NAME 

1  CWBGRID 


MESH  NAME:  PHANTGRID 


PEGSUS  MESH  NO.=  3 
PHANTOM  MESH 


INDEX  RANGES 

:  JMAX  = 

40  KMAX 

=  10 

LMAX  =  30 

TRANSLATION: 

ROTATION: 

XO 

XR 

ALPHA  = 

0 . 0000E+00 

0 . 0000E+00 

0 . OOOOE+OO 

Y0 

YR 

BETA  = 

0 . OOOOE+OO 

0. 0000E+00 

0. 0000E+00 

Z0  = 

ZR  = 

GAM  = 

0. 0000E+00 
0 . 0000E+00 
0. 0000E+00 

SCALING:  SCALE  =  0.1000E+01 


MAX  X  =  0 . 1000E+31 

MAX  Y  =  0 . 1000E+31 

MAX  Z  =  0 . 1000E+3 1 

MAX  J  =  43 

MAX  K  =  36 

MAX  L  =  34 


=  1.0000 
=  0.3000 

=  -0.1000 
=  0 . 1000E-02 
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ROTATION  SEQUENCE 


NO.  OF  FRINGE  POINTS 
QUALITY ( 1 ) 

QUALITY ( 2 ) 

QUALITY ( 3 ) 

EPS 

MIN/MAX  OF  MESH 
MIN  X  =  0 . 1100E+01 

MIN  Y  =  -0 . 1000E+00 

MIN  Z  =  0 . OOOOE+OO 


=  1.0000 
=  1.0000 
=  1.0000 
=  0 . 1000E-02 


MAX  X  =  0 . 1526E+01 
MAX  Y  =  0 . 1000E+00 
MAX  Z  =  0.3400E+00 


NO  INCLUDE  RANGES  SPECIFIED 


LINKS 

NO .  NAME 


*****  BOUNDARY  PARAMETER  SUMMARY  ***** 


BOUNDARY  NAME:  PHANT  HOLE  BOUND 

BOUNDARY  TYPE  =  HOLE 
CLOSED  =  T 

IS  PART  OF  MESH  =  PHANTGRID 


MIN/MAX  OF  HOLE  BOUNDARY 


MIN 

X  = 

0 . 1100E+01 

MAX  X  = 

0.1526E+01 

MIN 

Y  = 

-0 . 1000E+00 

MAX  Y  = 

0 . 1000E+00 

MIN 

Z  = 

0 . 0000E+00 

MAX  Z  = 

0 . 3400E+00 

MAKE  HOLE  IN: 

NO .  NAME 
1  CWBGRID 

BOUNDARY  NAME:  FIN  OUTER  BOUND 

BOUNDARY  TYPE  =  OUTER 

CLOSED  =  F 

IS  PART  OF  MESH  =  CFINGRID 

SURFACE  NO.  =  1 

ISPARTOF  :  PHANT  HOLE  BOUND 
SURFACE  RANGES 

MIN  J  =  1  MAX  J  =  40 

MIN  K  =  1  MAX  K  =  10 

MIN  L  =  30  MAX  L  =  30 

NVOUT  =  +L 

CONNECTING  MESH  =  PHANTGRID 

SURFACE  NO.  =  2 

ISPARTOF  :  PHANT  HOLE  BOUND 
SURFACE  RANGES 

MIN  J  =  1  MAX  J  =  40 
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MIN  K  = 

1 

MAX  K  =  10 

MIN  L  = 

1 

MAX  L  =  1 

NVOUT  =  -L 

CONNECTING 

MESH 

=  PHANTGRID 

SURFACE  NO. 

= 

3 

ISPARTOF  : 

PHANT  HOLE  BOUND 

SURFACE  RANGES 

MIN  J  = 

1 

MAX  J  =  1 

MIN  K  = 

1 

MAX  K  =  10 

MIN  L  = 

1 

MAX  L  =  30 

NVOUT  =  -J 

CONNECTING 

MESH 

=  PHANTGRID 

SURFACE  NO. 

- 

4 

ISPARTOF  : 

PHANT  HOLE  BOUND 

SURFACE  RANGES 

MIN  J  = 

40 

MAX  J  =  40 

MIN  K  = 

1 

MAX  K  =  10 

MIN  L  = 

1 

MAX  L  =  30 

NVOUT  =  +J 

CONNECTING 

MESH 

=  PHANTGRID 

SURFACE  NO. 

= 

5 

ISPARTOF  : 

PHANT  HOLE  BOUND 

SURFACE  RANGES 

MIN  J  = 

1 

MAX  J  =  40 

MIN  K  = 

1 

MAX  K  =  1 

MIN  L  = 

1 

MAX  L  =  30 

NVOUT  =  -K 

CONNECTING 

MESH 

=  PHANTGRID 

SURFACE  NO. 

= 

6 

ISPARTOF  : 

PHANT  HOLE  BOUND 

SURFACE  RANGES 

MIN  J  = 

1 

MAX  J  =  40 

MIN  K  = 

10 

MAX  K  =  10 

MIN  L  = 

1 

MAX  L  =  30 

NVOUT  =  +K 

CONNECTING 

MESH 

=  PHANTGRID 

SURFACE  NO, 

= 

7 

ISPARTOF  : 

FIN  i 

OUTER  BOUND 

SURFACE  RANGES 

MIN  J  = 

1 

MAX  J  =  43 

MIN  K  = 

1 

MAX  K  =  36 

MIN  L  = 

33 

MAX  L  =  34 

NVOUT  = 

CONNECTING 

MESH 

=  CFINGRID 

SURFACE  NO 

— 

8 

ISPARTOF  : 

FIN  < 

OUTER  BOUND 

SURFACE  RANGES 

MIN  J  = 

1 

MAX  J  =  2 

MIN  K  = 

1 

MAX  K  =  36 

MIN  L  = 

2 

MAX  L  =  32 

NVOUT  = 
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CONNECTING  MESH  =  CFINGRID 


SURFACE  NO. 

,  = 

9 

ISPARTOF  : 

FIN 

OUTER  BOUND 

SURFACE  RANGES 

MIN  J  = 

42 

MAX  J  = 

43 

MIN  K  = 

1 

MAX  K  = 

36 

MIN  L  = 

2 

MAX  L  = 

32 

NVOUT  = 

CONNECTING 

MESH 

=  CFINGRID 

SURFACE  NO, 

= 

10 

ISPARTOF  : 

FIN 

OUTER  BOUND 

SURFACE  RANGES 

MIN  J  = 

3 

MAX  J  = 

41 

MIN  K  = 

1 

MAX  K  = 

2 

MIN  L  = 

2 

MAX  L  = 

32 

NVOUT  = 

CONNECTING 

MESH 

=  CFINGRID 

SURFACE  NO. 

_ 

11 

ISPARTOF  : 

FIN 

OUTER  BOUND 

SURFACE  RANGES 

MIN  J  = 

3 

MAX  J  = 

41 

MIN  K  = 

35 

MAX  K  = 

36 

MIN  L  = 

2 

MAX  L  = 

32 

NVOUT  = 

CONNECTING 

MESH 

=  CFINGRID 

*****  REGION  PARAMETER  SUMMARY  ***** 

REGION  NAME:  FIN  REGION 

REGION  TYPE  =  INTR 
IS  PART  OF  MESH  =  CFINGRID 

VOLUME  NO.  =  1 

ISPARTOF  :  FIN  REGION 


VOLUME  RANGES 


MIN 

J  = 

14 

MAX  J  = 

29 

MIN 

K  = 

18 

MAX  K  = 

19 

MIN 

L  = 

2 

MAX  L  = 

20 

INTERPOLATING  BOUNDARY  POINTS  FOR 
MESH:  CWBGRID 
LINK:  CFINGRID 

QUALITY=  1.00000  STENCILS  FOUND: 

NO  ORPHAN  POINTS  LEFT  IN  MESH  -  CWBGRID 

INTERPOLATING  BOUNDARY  POINTS  FOR 
MESH:  CFINGRID 
LINK:  CWBGRID 

QUALITY=  1.00000  STENCILS  FOUND: 

LINK:  CWBGRID 
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1712 


12210 


QUALITY=  0.90000 
LINK:  CWBGRID 

STENCILS 

FOUND: 

76 

QUALITY =  0.80000 
LINK:  CWBGRID 

STENCILS 

FOUND: 

76 

QUALITY=  0.70000 
LINK:  CWBGRID 

STENCILS 

FOUND: 

20 

QUALITY=  0.60000 
LINK:  CWBGRID 

STENCILS 

FOUND: 

12 

QUALITY=  0.50000 

STENCILS 

FOUND: 

2 

NO  ORPHAN  POINTS  LEFT 

IN  MESH 

-  CFINGRID 

INTERPOLATING 

BOUNDARY 

POINTS 

FOR 

MESH:  PHANTGRID 

*****  QUALITY  SUMMARY 

***** 

MINIMUM 

NO.  OF  INTERPOLATED 

QUALITY 

BOUNDARY 

POINTS 

FOUND 

=  1.0 

13922 

>  =  0.9 

76 

00 

o 

it 

A 

76 

O 

o 

II 

A 

20 

O 

II 

A 

12 

>=0.5 

2 

o 

II 

A 

0 

co 

o 

II 

A 

0 

V 

II 

o 

ES-) 

0 

V 

II 

o 

I—4 

0 

V 

II 

o 

o 

0 

MESH:  CWBGRID 

NUMBER  OF  BOUNDARY  POINTS:  1712 

NUMBER  OF  INTERPOLATION  STENCILS:  12396 

MESH:  CFINGRID 

NUMBER  OF  BOUNDARY  POINTS:  12396 

NUMBER  OF  INTERPOLATION  STENCILS:  1712 


XMER3D  PARAMETERS 


MINIMUM 

DIMENSIONS 

REQUIRED 

FOR 

THE  FLOW 

SOLVER 

IWRK 

=  385128 

IGRD 

- 

179725 

JDIM 

=  79 

KDIM 

= 

36 

LDIM 

65 

IIDIM 

=  12396 

IBDIM 

= 

12396 

IQBDIM  = 

14108 

XLIM3D  PARAMETERS 

MINIMUM 

DIMENSIONS 

REQUIRED 

FOR 

THE  FLOW 

SOLVER 

IWRK 

=  61992 

IGRD 

= 

179725 

LARG 

=  79 

NORPH 

= 

0 

IIDIM 

=  12396 

IBDIM 

= 

12396 

ISDIM  = 

898625 
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B2  Diagnostic  Maps 

DIAGNOSTIC  MAPS 

LEGEND 

BOUNDARY  POINTS  UPDATED  BY  : 

A  -  CWBGRID 
B  -  CFINGRID 

INTERPOLATION  STENCILS  UPDATING  POINTS  IN  : 
a  -  CWBGRID 
b  -  CFINGRID 

.  -  FIELD  POINT 

#  -  HOLE  POINT 

?  -  ORPHANED  BOUNDARY  POINT 

*  -  INTERPOLATION  STENCIL  UPDATING  MORE  THAN  ONE 

INTERPOLATION  BOUNDARY  POINT 

Here  in  the  following  are  samples  of  the  diagnostic  maps  for  the  two  meshes 

MAP  FOR  CWBGRID 

PLANE  L  =  31 

J12345678 

12345678901234567890123456789012345678901234567890123456789012345678901234567890 

K 

1  . * . ################# . 

2  . *b _ BB#################BB.  .** . 

3  . ** _ BB#################BB.  .** . 

4  . ** _ BB#################BB.  .** . 

5  . * _ BB#################BB.  .** . 

6  . * . BB################BB.  .** . 

7  . ** _ BB################BB  .  .** . 

8  . b* _ BB################BB.  .  *  * . 

9  . **. . . BB################BB . .** . 

10  . ** _ BB###############BB.  .** . 

11  . *b.  .  .  BB###############BB  .  .** . 

12  . ** _ BB##  ##########  ##BB  .  .** . 

13  . **.  .  .BB##############BB.  .** . 

14  . ** _ BB#############BB.  .** . 

15  . **.  .  .BB#############BB.  .** . 

16  . **.  .  .BB############BB.  .** . 

17  . . . **.  .  .BB###########BB.  .** . 

18  . **.  .  .BB##########BB.  .** . 

19  . **.  .  .BB#########BB.  .** . 

20  . *  * .  .  . BBBBBBB###BB. . ** . 

21  . ********  *bBBBBB .  .  .  *  * . 

22  . ************B. .  .  .  ** . 

23  . ********* . 

24  . 

25  . 

26  . 

27  . 

28  . 

29  . 

30  . 

31  . 

32  . 

33  . 

34  . 

35  . 
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MAP  FOR  CWBGRID 


PLANE  L  =  32 

J12345678 

12345678901234567890123456789012345678901234567890123456789012345678901234567890 

K 

1  . * . ################# . 

2  . ** _ BB#################BB .  .** . 

3  . ** _ BB#################BB.  .** . 

4  . ** _ BB#################BB.  .** . 

5  . * . BB################BB .  .** . 

6  . * . BB################BB .  .** . 

7  ** _ BB################BB .  .** . 

8  . ** _ BB################BB .  .** . 

9  . *b _ BB###############BB  .  .** . 

10  . ** _ BB###############BB.  .** . 

11  . * _ BB###############BB.  .** . 

12  . ** _ BB##############BB.  .** . 

13  . * _ BB##############BB.  .** . 

14  . ** _ BB#############BB.  .** . 

15  . **.  .  .BB#############BB.  .** . 

16  . **. . .BB############BB. .** . 

17  . **.  .  .BB###########BB.  .** . 

18  . **.  .  .BB##########BB.  .** . 

19  . *.  .  .  .  BBBBBBB####BB  .  .** . 

20  . *  *  .  .  .  BBBBBBBB##BB .  .  *  * . 

21  . **********  BBBBB .  .  .  *  * . 

22  . ************B - ** . 

23  . ******** . 

24  . 

25  . 

26  . 

27  . 

28  . 

29  . 

30  . 

31  . 

32  . 

33  . 

34  . 

35  . 
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MAP  FOR  CWBGRID 


PLANE  L  =  33 

J12345678 

12345678901234567890123456789012345678901234567890123456789012345678901234567890 

K 

1  . * . ################# . 

2  *  *  * . ** _ BB#################BB.  .** . 

3  . ** _ BB#################BB.  .** . 

4  . ** _ BB#################BB  .  .  ** . 

5  . * . BB################BB  .  .** . 

6  . * . BB################BB  .  .** . 

7  . ** _ BB################BB.  .** . 

8  . ** _ BB################BB.  .** . 

9  . *b _ BB###############BB .  .** . 

10  . ** _ BB###############BB.  .** . 

11  . * _ BB###############BB  .  .**.  . . . 

12  . ** _ BB##############BB .  .** . 

13  . * _ BB##############BB  .  .** . 

14  . ** _ BB# ########### #BB .  .** . 

15  . **.  .  ,BB#############BB.  .** . 

16  . **.  .  .BB############BB.  .** . 

17  . **.  .  .  BB###########BB  .  .** . 

18  . **.  .  .BB##########BB.  .** . 

19  . *  .  .  .  .BBBBBBB####BB.  .  ** . 

20  . **.  .  .  BBBBBBBB##BB  .  .** . 

21  . ******* ***BBBBB.  ..  ** . 

22  . ************8..,.** . 

23  . ******** . 

24  . 

25  . 

26  . 

27  . 

28  . 

29  . 

30  . 

31  . 

32  . 

33  . 

34  . 

35  . 
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MAP  FOR  CWBGRID 


PLANE  L  =  34 

J12345678 

12345678901234567890123456789012345678901234567890123456789012345678901234567890 

K 

1  . * . ################# . 

2  *b _ BB#################BB.  .** . 

3  ** _ BB#################BB.  .** . 

4  ”  *  * . ** _ BB### # ########### ##BB  .  .** . 

5  . * . BB################BB.  .** . 

6  . * . BB################BB.  .** . 

7  ** _ BB################BB.  .** . 

8  . b* _ BB#### ########### #BB .  .  ** . 

9  . ** _ BB###############BB  .  .** . 

10  . ** _ BB# # ############ #BB .  .** . 

11  . *b.  . . BB###############BB.  .** . 

12  . ** _ BB#### ###### ####BB .  .** . 

13  . **.  .  .BB##############BB.  .** . 

14  . ** _ BB#############BB  .  .** . 

15  . **.  .  .BB#############BB.  .** . 

IS  . **.  .  .BB############BB.  .** . 

17  . **..  .  BB###########BB.  .** . 

18  !! . **.  .  .BB##########BB.  .** . 

19  . *  *  .  .  .  BBBBBBB#  #  #  #BB  .  .  *  * . 

20  . *  *  .  .  .  BBBBBBBB##BB  .  .  *  * . 

21  . ******** *bBBBBB  .  .  .  *  * . 

22  . ************B.  .  .  .  ** . 

23  . ********* . 

24  . 

25  . 

26  . 

27  . 

28  . 

29  . 

30  . 

31  . 

32  . 

33  . 

34  . 

35  . 
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MAP  FOR  CWBGRID 


PLANE  L  =  35 

J12345678 

12345678901234567890123456789012345678901234567890123456789012345678901234567890 

K 

1  . * . ################# . 

2  . * . BB#################BB.  .** . 

3  . ** _ BB#################BB  .  .** . 

4  . b* _ BB#################BB.  .** . 

5  . * _ BB# ########### #####BB.  .** . 

6  . ** _ BB################BB .  .** . 

7  . ** _ BB################BB.  .** . 

8  . *b. . - BB################BB . .** . 

9  . **.  .  . BB################BB .  .** . 

10  . b* _ BB###############BB  .  .** . 

11  . **.  .  .  BB###############BB  .  .** . 

12  . b* _ BB##############BB .  .** . 

13  . **.  .  .BB##############BB.  .** . 

14  . * _ BB#############BB.  .** . 

15  . **.  .  .BB#############BB.  .** . 

16  . **.  .  .BB############BB.  .** . 

17  . **.  .  .BB###########BB.  .** . 

18  . **.  .  .BB##########BB.  .** . 

19  . **.  .  .BB#########BB.  .** . 

20  . *  .  .  .  .  BBBBBBB###BB  .  .  ** . 

21  . *  *  .  .  .  BBBBBBBBBB  .  .  .  *  * . 

22  !!! . *********  *bBB  . 

23  . k*********** . 

24  . 

25  . 

26  . 

27  . 

28  . 

29  . 

30  . 

31  . 

32  . 

33  . 

34  . 

35  . 
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MAP  FOR  CFINGRID 


PLANE  L  =  2 


J  1  2  3  4  5 

12345678901234567890123456789012345678901234567890 


1  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

2  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


3  AA . 

4  AA . 

5  AA . a . a . 

6  AA . aaaaaaaa .  a .  aaa .  a . 

7  AA . a.aa.  .a.a . 

8  AA . aaa . a*.aa.a . 

9  AA . 

10  AA. . . . aa . 

11  AA . aa.  .  .  . 

12  AA. . . . aa . 

13  AA . 

14  AA . a.a . aa - 

15  AA . 

1 6  AA . a.a . aa .  .  .  . 

17  AA . 

18  AA . a.a...################ . aa - 

19  AA . ################ . 

20  AA . a.a . aa .  .  .  . 

21  AA . 

22  AA . a.a . aa.  .  .  . 

23  AA . 

24  AA. . . . aa . 

25  AA . aa - 

26  AA. . . .aa . 

27  AA . 

28  AA . aaa . a*  .aa.  a . 

2  9  AA . a .  aa  ..a.a . 

30  AA . aaaaaaaa .  a .  aaa .  a . 

31  AA . a . a . 

32  AA . 

33  AA . 

34  AA . 


35  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

36  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
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MAP  FOR  CFINGRID 


PLANE  L  = 


J  1  2  3  4  5 

1234567890X234567890123456789012345678901234567890 


1  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

2  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

3  AA . AA 

4  AA . AA 

5  AA . a* . AA 

6  AA . aaaaaa . aa . AA 

7  AA . a .  aa.  .  a.  a . AA 

8  AA . aaaaaa.a.aaa.aa.a . AA 

9  AA . a..a..a.a.aa . AA 

10  AA.  .  .  .a . AA 

11  AA . aa.  .  .  .AA 

12  AA - aa . AA 

13  AA . AA 

14  AA . aa . aa.  .  .  . AA 

15  AA . AA 

16  AA . aa . aa.  .  .  . AA 

17  AA . AA 

18  AA . aa _ ################ . aa - AA 

19  AA . ################ . AA 

20  AA . aa . aa.  .  .  .AA 

21  AA . AA 

22  AA . aa . aa.  .  .  .AA 

23  AA . AA 

24  AA.  .  .  .aa . AA 

25  AA . aa.  .  .  . AA 

26  AA _ a . AA 

27  AA . a..a..a.a.aa . AA 

28  AA . aaaaaa.a.aaa.aa.a . AA 

2  9  AA . a .  aa .  .  a .  a . AA 

30  AA . aaaaaa . aa . AA 

31  AA . a* . AA 

32  AA . AA 

33  AA . AA 

34  AA . AA 

35  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

36  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
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MAP  FOR  CFINGRID 


PLANE  L  = 


J  1  2  3  4  5 

12345678901234567890123456789012345678901234567890 


1  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

2  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

3  AA . AA 

4  AA . AA 

5  AA . a.aa . AA 

6  AA . a. aaaaaaaa . a . aa . AA 

7  AA . a. a. .a . a . AA 

8  AA . a.  aaaaaaa .  a .  . aa . aa . a . AA 

9  AA . aa.a.a.  .a . AA 

10  AA . AA 

11  AA . aa . aa - AA 

12  AA . AA 

13  AA . a . AA 

14  AA . aa . a - AA 

15  AA . AA 

16  AA . aa . aa.  .  .  . AA 

17  AA . AA 

18  AA . aa _ ################ . aa - AA 

19  AA . ################ . AA 

20  AA . aa . aa...  .AA 

21  AA . AA 

22  AA . aa . a.  .  .  .AA 

23  AA . a . AA 

24  AA . AA 

25  AA . aa . aa - AA 

26  AA . AA 

27  AA . aa.a.a.  .a . AA 

28  AA . a.  aaaaaaa. a.  .aa.aa.  a . AA 

29  AA . a. a.  .a . a . AA 

30  AA . a. aaaaaaaa . a . aa . AA 

31  AA . a.aa . AA 

32  AA . AA 

33  AA . AA 

34  AA . AA 

35  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
3  6  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
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MAP  FOR  CFINGRID 


PLANE  L  =  33 


J  1  2  3  4  5 

12345678901234567890123456789012345678901234567890 

K 

1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

2  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

3  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

4  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

5  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

6  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

7  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

8  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

9  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

10  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

11  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

12  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

13  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

14  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

15  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

16  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

17  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

18  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

19  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

20  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

21  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

22  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

23  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

24  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

25  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

26  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

27  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

28  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

29  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

30  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

31  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

32  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

33  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

34  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

35  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

36  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
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MAP  FOR  CFINGRID 


PLANE  L  =  34 


J  1  2  3  4  5 

12345678901234567890123456789012345678901234567890 
K 

1  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

2  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

3  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

4  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

5  AAAlAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

6  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAuAAAAAAAAAA 

7  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

8  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

9  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAuAA 

10  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

11  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

12  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

13  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

14  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

15  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

16  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

17  AAAAAuAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

18  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

19  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

20  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

21  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

22  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

23  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

24  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

25  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

26  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

27  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

28  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

29  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
3  0  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

31  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAlAAAAAAAAAAAAA 

32  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

33  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

34  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
3  5  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
36  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
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Appendix  C  -  Numerical  Application  of  the 
Euler  Boundary  Conditions 


The  classical  Euler  boundary  conditions  have  been  known  for  decades  as  “there  is  no 
penetration  of  the  fluid  medium  to  the  solid  boundary”,  or  alternatively  the  normal  component  of 
the  relative  fluid  velocity  is  zero  on  a  solid  boundary.  Although  the  statement  of  the  Euler 
boundary  condition  seems  simple,  its  application  in  computational  fluid  dynamics  is  not.  In  the 
rest  of  this  Appendix,  the  numerical  application  of  the  Euler  boundary  conditions  is  presented  in 
detail.  The  methodology  will  be  discussed  for  a  general  case  when  the  grids  are  allowed  to  move, 
and  the  static  case  can  be  easily  derived  by  assigning  zero  value  for  the  grid  velocity. 

The  following  assumptions  have  been  made  in  the  next  derivation: 

First-order  approximation  for  the  tangents  at  the  grid  points  on  the  body  surface. 

The  magnitude  of  the  velocity  at  a  grid  point  on  the  body  surface  is  interpolated  from  the  next 
point  to  the  surface  in  the  direction  of  the  outward  normal. 


Figure  C-l.  Schematic  drawing  for  the  body  surface. 

Let  Vj  and  v2  denote  the  two  vectors  that  representing  the  first-order  approximation  of  the 
tangents  to  the  surface  at  the  grid  point  (see  figure  C-l). 
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V,  =  (x O'  + 1,  j,  k )  -  x(i  - 1,  j,  k))i  +  ( y(i  + 1,  j,k)  -  y{i  - 1  ,j,  k))j  +  (z(i  +  l,j,k)~  z(i  - 1,  j,  k))k 

=  vlJ  +  \J  +  \k 

v2  =  (x (i,  j  +  U)  -  x(i,  j  - 1,  k))i  +  (y(i,  j  +  U)  -  y(i,  j  -  l kj)j  +  (z(i,  j  +  U)  -  z(i,  j  - 1,  k))k 
=  v2J  +  v2J  +  v2k 

The  outward  normal  to  the  surface  N  is  determined  as  follows: 

N  =  v,  x  v2 

=  (viA  -  v,,v2y )/  +  (v1;v2x  - vljV2i );  +  (vIxv2y  -  vIyv2x  )k 
=  Nxi  +  Nyj  +  Nzk 

Since  v,  and  v2  are  not  orthogonal  in  general,  let  tx  and  n  denote  the  unit  vectors  in  the 
direction  of  v,  and  N  respectively,  then  use  tx  and  n  to  find  the  third  orthogonal  unit  vector  t2 . 
t2  =  n  x  F, 

Where: 


?  *  f 

=  tx  i  +tx  j  +  tx  k 

_  N 


=  nxi  +  nyj  +  n.k 

/V  ^  ^ 

t2  =  (Vi,  _  nzf i, ) *  +  ("A  ~  »  A )./  +  («A  -  «  A  )* 

A  A  A 

—  f2  i  +  f2  j  +  f2  ^ 
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Assume  the  body  is  moving  with  a  velocity  U b ,  where: 


Ub  -  uj  +  vb]  +  wbk 

Then  the  body  velocity  in  the  normal  direction  is  given  by: 
ub„  =  Ub-n 

=  ubnx+vbny+wbnz 

Projecting  the  fluid  velocity  U  on  the  body  surface  along  the  three  orthogonal  directions  n,tx, 
and  t2  yields: 

un  =U  -n 


=  unx  +vny  +  wnz 
u.  =  V  ■  t, 

M  1 

-  u  t,  +  vr, 

Ay  z 

uh  =  V -t2 


-ut2  +Vt2,  +Wt2 

Now,  apply  the  Euler  boundary  condition: 

un  =  ub„  (C-1) 

Interpolate  the  value  of  the  velocity  components  in  the  other  two  directions  from  that  of  the  next 
layer  in  the  direction  of  the  outward  normal: 


(C-2) 

(C-3) 


The  ±  sign  indicates  that,  the  interpolation  will  be  done  from  either  the  next  or  previous  layer 
according  to  the  direction  of  the  outward  normal.  For  the  directions  posted  in  figure  C-l,  the 
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outward  normal  is  running  in  the  direction  of  increasing  K,  so  the  interpolation  will  be  done  from 
the  surface  ( K+l )  and  vice-versa. 

Equations  (C-l  to  C-3)  could  be  written  in  matrix  form  as  follows: 


n 

n„ 

(u ' 

X 

y 

z 

f. 

h 

h 

V 

lx 

h 

t2 

*2Z 

w 

V  / 

ubnx+v„ny+wbnt 

t]  U±  +fj  V±  +?!  W± 
^2  W+  +  £2  V+  ~^2 


(C-4) 


Equation  (C-4)  is  a  system  of  three  linear  determined  algebraic  equations,  which  can  be  easily 
solved  by  Kramer’s  Rule.  Finally,  the  solution  is  given  by: 


u  ~  Ub„  ^2Z 

V  =  nMj2l 
w  =  nAu,Jly 


- \ *2, ) - »>•  («,it hz-\uH±)+nz K \  - \ uhi ) 

~^lzUt2A~Ubn^lj2z  ~  h^2x  )  "*■  Uz  (Ut2i *1,  ~*2xUt,A 

-  t2y  Ut{i)~ny  U‘z±  ~  U,I±  *2  A +  Ub.^lx  ^2,  —  *1,  ^2X  ) 
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Appendix  D  —  Numerical  Computation  of  the 
Aerodynamic  Forces  and  Moments 

Given  an  arbitrary  body  immersed  in  a  fluid  medium  (figure  D-l),  the  elementary  force 
dF  exerted  on  an  element  of  the  surface  area  do  is  given  by: 

dF  =  (n-f)do 

Where: 

ti  Denotes  the  outward  normal  to  the  surface, 
f  Denotes  the  stress  tensor. 


Figure  D-l.  Schematic  drawing  shows  the  nomenclature  used. 

Let  the  equation  of  the  surface  be  given  by: 

f(x,y,z)  =  T](x,y,z)  =  c  cel 

Then  by  definition  the  normal  vector  h  is  given  by: 


Vj  +  Vyj  +  V^ 

Jnf+vf+vf 


dFx  =  -P  rfx  do 
dFy  =  -P  ffy  do 
dFz  =  -P  rfz  do 
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Assume  dSx  and  dS2  are  two  arbitrary  orthogonal  curvilinear  elements  of  arc  length  that  bound 
the  element  of  area  do  .  For  simplicity,  let  us  pick  dSx  along  the  line  £  =  const,  and  dS2  along 
the  line  £  =  const. 

In  general,  |r£S|  =  -Jdx~  +dy~  +  dz~ 

Also  dx,  dy,  and  dz  could  be  written  as: 

dx  =  x~  dE,  +  x^dl]  +  a-,  dE, 

dy  =  y(d^  +  yndt]  +  y(d£ 
dz  -  z^dE,  +  zndri  +  z^dC, 

Substitute  by  dE,  =0on  dS{  and  dE,  =0on  dS2.  Since  dSx  and  dS2  are  bounding  the  element 
of  area  do  ,  then  drj  is  also  zero  on  dSx  and  dS2 . 

=>  \dSx  |  =  J(x(  d£  J  +  (yc  d£  f  +  (z(d£  f 

c?52|  =  -J(x^dZj  +  (y^d^f  +(z€dt;f 

By  definition,  do  is  given  by: 

do  =  dSx  x  dS2 

=  jjc/S,  | 2 1 dS2 1 2  -  (dSx  ■  dS2  J  ^ 

Factoring  out  the  previous  equation,  and  inserting  the  expressions  for  dSx  and  dS2 . 

=>  d<7  =  y |  —  Xg  y  ^  +  (x^  z ^  —  y^z^f  d£d£  (D-2) 

Recall  the  following  metric  relations: 

nx  =  J{y(zi-y(zc) 

1,=J(xezc  ~xczt) 
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riz=j(x(y^-x^yc) 

Where  J  is  the  Jacobian  of  the  coordinate  transformation  at  the  grid  point  (£,?7,  £). 

Insert  the  metric  relations  into  (D-2),  finally,  do  is  given  by: 

do  =  yfil  +rfy  +Vl  (°-3) 

Inserting  (D-3)  into  (D-l)  yields: 

dF  =  -P^dtdt 
J 

=>  dFx=-P^d%dC 

dFy=-P^d£dC 

dF  =  -P—  dE,  dC, 
z  J 

The  differential  force  dF  exerts  a  moment  dM  about  the  center  of  gravity,  which  is  given  by: 

dM  =  rt  x  dF  (D-4) 

Where: 

r  Denotes  the  displacement  vector  from  the  c.g.  to  the  surface  area  do. 

r,  =(x-xcg)T  +  (y-ycg)j  +  (z-zcg)k 
Inserting  the  expression  for  ri  in  (D-4)  and  applying  the  cross  product  yields: 

dMx=(y-  ycg  )dFz-(z~  zcg  )dFy 
dM  y  ={z~  zcg  )  dFx  -(x-  xcg  )dFz 
dMz  =  (x  -  xcg )  dFy  -  (y  -  ycg  )dFx 
Inserting  the  expressions  for  dF x ,  dFy ,  and  dF.  yields: 


D-4 


dMx  =  -j [(y  -  ycg)t]z  - (z  -  zcg  )r]y \d^ d£ 
dMy  =-j[(z-  zcg)vx -O- *cs )VZ ]d£ dC 
dMz  =-j[(x-  xcg)r\y  -(y-  ycg )r]x ]d% d £ 

Integrating  over  the  whole  surface  yields  the  total  forces  and  moments  components,  which  are 
given  respectively  as  follows: 


J  P^dSdt 

Fz=-J 

]  P^didi 

F<=-J 

J  P'tfdtdt 

M,  =  -J. 

f  -ti-z„)n,]d£dt; 

M,=-J 

\  y[(z-z„)tj, -(x-Xt,K]d£d( 

M„=-J 

j  y  [(*  -  .v  ,  )?J,  -  (y  -  m  ,  \t!i 

Since  the  body-fixed  x-axis  is  only  translating  (not  rotating),  then  the  rolling  moment  about  the 
body-fixed  x-axis  (Lb)  in  equation  (3-3)  equals  Mx. 
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